Skip to content

Commit 2551fde

Browse files
DanielDoehringranocha
DanielDoehring
and
ranocha
authored
Type assertions indicators (#2643)
* fmt * MArray indicators * fmt * typo * 2D statically sized, above standard array * Update src/solvers/dgsem_tree/dg_1d.jl * Update src/solvers/dgsem_tree/indicators_3d.jl * Apply suggestions from code review * adapt * Update docs/src/conventions.md Co-authored-by: Hendrik Ranocha <[email protected]> * mvec * type assertions p4est --------- Co-authored-by: Hendrik Ranocha <[email protected]>
1 parent 8855d55 commit 2551fde

File tree

9 files changed

+112
-133
lines changed

9 files changed

+112
-133
lines changed

‎docs/src/conventions.md‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -125,6 +125,16 @@ based on the following rules.
125125
be used when necessary, e.g., to avoid additional overhead when interfacing
126126
with external C libraries such as HDF5, MPI, or visualization.
127127

128+
### Usage of statically sized arrays
129+
130+
Trixi.jl employs statically sized vectors and arrays from the
131+
[StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) package, i.e.,
132+
`SVector, SMatrix` for the immutable versions and `MVector, MMatrix` for the mutable variants.
133+
These types are preferred for "small" (StaticArrays.jl [recommends as a rule of thumb about 100 entries](https://juliaarrays.github.io/StaticArrays.jl/stable/#When-Static-Arrays-may-be-useful)) arrays with a size known at compile time.
134+
Inspired by this recommendation, Trixi.jl uses the static size `SVector` most notably for numerical fluxes and initial/boundary conditions.
135+
The mutable statically sized `MArray` type is used for **arrays up to two dimensions** which show up in mortars, indicators, and certain volume integrals.
136+
For arrays of higher dimensions (mostly three and four) we use standard Julia `Array` types.
137+
128138
## Numeric types and type stability
129139

130140
In Trixi.jl, we use generic programming to support custom data types to store the numerical simulation data, including standard floating point types and automatic differentiation types.

‎src/solvers/dgmulti/shock_capturing.jl‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -35,12 +35,13 @@ end
3535
# this method is used when the indicator is constructed as for shock-capturing volume integrals
3636
function create_cache(::Type{IndicatorHennemannGassner}, equations::AbstractEquations,
3737
basis::RefElemData{NDIMS}) where {NDIMS}
38-
alpha = Vector{real(basis)}()
38+
uEltype = real(basis)
39+
alpha = Vector{uEltype}()
3940
alpha_tmp = similar(alpha)
4041

41-
A = Vector{real(basis)}
42-
indicator_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.maxthreadid()]
43-
modal_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.maxthreadid()]
42+
MVec = MVector{nnodes(basis), uEltype}
43+
indicator_threaded = MVec[MVec(undef) for _ in 1:Threads.maxthreadid()]
44+
modal_threaded = MVec[MVec(undef) for _ in 1:Threads.maxthreadid()]
4445

4546
# initialize inverse Vandermonde matrices at Gauss-Legendre nodes
4647
(; N) = basis

‎src/solvers/dgsem_p4est/dg_3d.jl‎

Lines changed: 15 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -10,21 +10,21 @@
1010
function create_cache(mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations,
1111
mortar_l2::LobattoLegendreMortarL2, uEltype)
1212
# TODO: Taal compare performance of different types
13-
fstar_primary_threaded = [Array{uEltype, 4}(undef, nvariables(equations),
14-
nnodes(mortar_l2),
15-
nnodes(mortar_l2), 4)
16-
for _ in 1:Threads.maxthreadid()] |> VecOfArrays
17-
fstar_secondary_threaded = [Array{uEltype, 4}(undef, nvariables(equations),
18-
nnodes(mortar_l2),
19-
nnodes(mortar_l2), 4)
20-
for _ in 1:Threads.maxthreadid()] |> VecOfArrays
21-
22-
fstar_tmp_threaded = [Array{uEltype, 3}(undef, nvariables(equations),
23-
nnodes(mortar_l2), nnodes(mortar_l2))
24-
for _ in 1:Threads.maxthreadid()] |> VecOfArrays
25-
u_threaded = [Array{uEltype, 3}(undef, nvariables(equations), nnodes(mortar_l2),
26-
nnodes(mortar_l2))
27-
for _ in 1:Threads.maxthreadid()] |> VecOfArrays
13+
A4d = Array{uEltype, 4}
14+
fstar_primary_threaded = A4d[A4d(undef, nvariables(equations),
15+
nnodes(mortar_l2), nnodes(mortar_l2), 4)
16+
for _ in 1:Threads.maxthreadid()] |> VecOfArrays
17+
fstar_secondary_threaded = A4d[A4d(undef, nvariables(equations),
18+
nnodes(mortar_l2), nnodes(mortar_l2), 4)
19+
for _ in 1:Threads.maxthreadid()] |> VecOfArrays
20+
21+
A3d = Array{uEltype, 3}
22+
fstar_tmp_threaded = A3d[A3d(undef, nvariables(equations),
23+
nnodes(mortar_l2), nnodes(mortar_l2))
24+
for _ in 1:Threads.maxthreadid()] |> VecOfArrays
25+
u_threaded = A3d[A3d(undef, nvariables(equations),
26+
nnodes(mortar_l2), nnodes(mortar_l2))
27+
for _ in 1:Threads.maxthreadid()] |> VecOfArrays
2828

2929
(; fstar_primary_threaded, fstar_secondary_threaded, fstar_tmp_threaded, u_threaded)
3030
end

‎src/solvers/dgsem_tree/dg_1d.jl‎

Lines changed: 4 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -40,13 +40,10 @@ function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations,
4040
VolumeIntegralFluxDifferencing(volume_integral.volume_flux_dg),
4141
dg, uEltype)
4242

43-
A2d = Array{uEltype, 2}
44-
fstar1_L_threaded = A2d[A2d(undef, nvariables(equations),
45-
nnodes(dg) + 1)
46-
for _ in 1:Threads.maxthreadid()]
47-
fstar1_R_threaded = A2d[A2d(undef, nvariables(equations),
48-
nnodes(dg) + 1)
49-
for _ in 1:Threads.maxthreadid()]
43+
MA2d = MArray{Tuple{nvariables(equations), nnodes(dg) + 1},
44+
uEltype, 2, nvariables(equations) * (nnodes(dg) + 1)}
45+
fstar1_L_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
46+
fstar1_R_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
5047

5148
return (; cache..., fstar1_L_threaded, fstar1_R_threaded)
5249
end

‎src/solvers/dgsem_tree/dg_2d.jl‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -96,10 +96,6 @@ function create_cache(mesh::TreeMesh{2}, equations,
9696
fstar_secondary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
9797
fstar_secondary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
9898

99-
# A2d = Array{uEltype, 2}
100-
# fstar_upper_threaded = [A2d(undef, nvariables(equations), nnodes(mortar_l2)) for _ in 1:Threads.maxthreadid()]
101-
# fstar_lower_threaded = [A2d(undef, nvariables(equations), nnodes(mortar_l2)) for _ in 1:Threads.maxthreadid()]
102-
10399
cache = (; fstar_primary_upper_threaded, fstar_primary_lower_threaded,
104100
fstar_secondary_upper_threaded, fstar_secondary_lower_threaded)
105101

‎src/solvers/dgsem_tree/dg_3d.jl‎

Lines changed: 18 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -45,23 +45,23 @@ function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3},
4545
volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG,
4646
uEltype)
4747
A4d = Array{uEltype, 4}
48-
fstar1_L_threaded = A4d[A4d(undef, nvariables(equations), nnodes(dg) + 1,
49-
nnodes(dg), nnodes(dg))
48+
fstar1_L_threaded = A4d[A4d(undef, nvariables(equations),
49+
nnodes(dg) + 1, nnodes(dg), nnodes(dg))
5050
for _ in 1:Threads.maxthreadid()]
51-
fstar1_R_threaded = A4d[A4d(undef, nvariables(equations), nnodes(dg) + 1,
52-
nnodes(dg), nnodes(dg))
51+
fstar1_R_threaded = A4d[A4d(undef, nvariables(equations),
52+
nnodes(dg) + 1, nnodes(dg), nnodes(dg))
5353
for _ in 1:Threads.maxthreadid()]
54-
fstar2_L_threaded = A4d[A4d(undef, nvariables(equations), nnodes(dg),
55-
nnodes(dg) + 1, nnodes(dg))
54+
fstar2_L_threaded = A4d[A4d(undef, nvariables(equations),
55+
nnodes(dg), nnodes(dg) + 1, nnodes(dg))
5656
for _ in 1:Threads.maxthreadid()]
57-
fstar2_R_threaded = A4d[A4d(undef, nvariables(equations), nnodes(dg),
58-
nnodes(dg) + 1, nnodes(dg))
57+
fstar2_R_threaded = A4d[A4d(undef, nvariables(equations),
58+
nnodes(dg), nnodes(dg) + 1, nnodes(dg))
5959
for _ in 1:Threads.maxthreadid()]
60-
fstar3_L_threaded = A4d[A4d(undef, nvariables(equations), nnodes(dg),
61-
nnodes(dg), nnodes(dg) + 1)
60+
fstar3_L_threaded = A4d[A4d(undef, nvariables(equations),
61+
nnodes(dg), nnodes(dg), nnodes(dg) + 1)
6262
for _ in 1:Threads.maxthreadid()]
63-
fstar3_R_threaded = A4d[A4d(undef, nvariables(equations), nnodes(dg),
64-
nnodes(dg), nnodes(dg) + 1)
63+
fstar3_R_threaded = A4d[A4d(undef, nvariables(equations),
64+
nnodes(dg), nnodes(dg), nnodes(dg) + 1)
6565
for _ in 1:Threads.maxthreadid()]
6666

6767
return (; fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded,
@@ -80,14 +80,16 @@ function create_cache(mesh::TreeMesh{3}, equations,
8080
nnodes(mortar_l2))
8181
for _ in 1:Threads.maxthreadid()]
8282
fstar_primary_upper_right_threaded = A3d[A3d(undef, nvariables(equations),
83-
nnodes(mortar_l2), nnodes(mortar_l2))
83+
nnodes(mortar_l2),
84+
nnodes(mortar_l2))
8485
for _ in 1:Threads.maxthreadid()]
8586
fstar_primary_lower_left_threaded = A3d[A3d(undef, nvariables(equations),
8687
nnodes(mortar_l2),
8788
nnodes(mortar_l2))
8889
for _ in 1:Threads.maxthreadid()]
8990
fstar_primary_lower_right_threaded = A3d[A3d(undef, nvariables(equations),
90-
nnodes(mortar_l2), nnodes(mortar_l2))
91+
nnodes(mortar_l2),
92+
nnodes(mortar_l2))
9193
for _ in 1:Threads.maxthreadid()]
9294
fstar_secondary_upper_left_threaded = A3d[A3d(undef, nvariables(equations),
9395
nnodes(mortar_l2),
@@ -105,7 +107,8 @@ function create_cache(mesh::TreeMesh{3}, equations,
105107
nnodes(mortar_l2),
106108
nnodes(mortar_l2))
107109
for _ in 1:Threads.maxthreadid()]
108-
fstar_tmp1_threaded = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2),
110+
fstar_tmp1_threaded = A3d[A3d(undef, nvariables(equations),
111+
nnodes(mortar_l2),
109112
nnodes(mortar_l2))
110113
for _ in 1:Threads.maxthreadid()]
111114

‎src/solvers/dgsem_tree/indicators_1d.jl‎

Lines changed: 18 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -8,19 +8,21 @@
88
# this method is used when the indicator is constructed as for shock-capturing volume integrals
99
function create_cache(::Type{IndicatorHennemannGassner},
1010
equations::AbstractEquations{1}, basis::LobattoLegendreBasis)
11-
alpha = Vector{real(basis)}()
11+
uEltype = real(basis)
12+
alpha = Vector{uEltype}()
1213
alpha_tmp = similar(alpha)
1314

14-
A = Array{real(basis), ndims(equations)}
15-
indicator_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.maxthreadid()]
16-
modal_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.maxthreadid()]
15+
MVec = MVector{nnodes(basis), uEltype}
16+
17+
indicator_threaded = MVec[MVec(undef) for _ in 1:Threads.maxthreadid()]
18+
modal_threaded = MVec[MVec(undef) for _ in 1:Threads.maxthreadid()]
1719

1820
return (; alpha, alpha_tmp, indicator_threaded, modal_threaded)
1921
end
2022

2123
# this method is used when the indicator is constructed as for AMR
22-
function create_cache(typ::Type{IndicatorHennemannGassner}, mesh,
23-
equations::AbstractEquations, dg::DGSEM, cache)
24+
function create_cache(typ::Type{IndicatorHennemannGassner},
25+
mesh, equations::AbstractEquations, dg::DGSEM, cache)
2426
create_cache(typ, equations, dg.basis)
2527
end
2628

@@ -111,20 +113,22 @@ function apply_smoothing!(mesh::TreeMesh{1}, alpha, alpha_tmp, dg, cache)
111113
end
112114

113115
# this method is used when the indicator is constructed as for shock-capturing volume integrals
114-
function create_cache(::Type{IndicatorLöhner}, equations::AbstractEquations{1},
115-
basis::LobattoLegendreBasis)
116-
alpha = Vector{real(basis)}()
116+
function create_cache(::Union{Type{IndicatorLöhner}, Type{IndicatorMax}},
117+
equations::AbstractEquations{1}, basis::LobattoLegendreBasis)
118+
uEltype = real(basis)
119+
alpha = Vector{uEltype}()
117120

118-
A = Array{real(basis), ndims(equations)}
119-
indicator_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.maxthreadid()]
121+
MVec = MVector{nnodes(basis), uEltype}
122+
123+
indicator_threaded = MVec[MVec(undef) for _ in 1:Threads.maxthreadid()]
120124

121125
return (; alpha, indicator_threaded)
122126
end
123127

124128
# this method is used when the indicator is constructed as for AMR
125-
function create_cache(typ::Type{IndicatorLöhner}, mesh, equations::AbstractEquations,
126-
dg::DGSEM, cache)
127-
create_cache(typ, equations, dg.basis)
129+
function create_cache(typ::Union{Type{IndicatorLöhner}, Type{IndicatorMax}},
130+
mesh, equations::AbstractEquations, dg::DGSEM, cache)
131+
return create_cache(typ, equations, dg.basis)
128132
end
129133

130134
function (löhner::IndicatorLöhner)(u::AbstractArray{<:Any, 3},
@@ -160,23 +164,6 @@ function (löhner::IndicatorLöhner)(u::AbstractArray{<:Any, 3},
160164
return alpha
161165
end
162166

163-
# this method is used when the indicator is constructed as for shock-capturing volume integrals
164-
function create_cache(::Type{IndicatorMax}, equations::AbstractEquations{1},
165-
basis::LobattoLegendreBasis)
166-
alpha = Vector{real(basis)}()
167-
168-
A = Array{real(basis), ndims(equations)}
169-
indicator_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.maxthreadid()]
170-
171-
return (; alpha, indicator_threaded)
172-
end
173-
174-
# this method is used when the indicator is constructed as for AMR
175-
function create_cache(typ::Type{IndicatorMax}, mesh, equations::AbstractEquations,
176-
dg::DGSEM, cache)
177-
cache = create_cache(typ, equations, dg.basis)
178-
end
179-
180167
function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 3},
181168
mesh, equations, dg::DGSEM, cache;
182169
kwargs...)

‎src/solvers/dgsem_tree/indicators_2d.jl‎

Lines changed: 16 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -8,16 +8,16 @@
88
# this method is used when the indicator is constructed as for shock-capturing volume integrals
99
function create_cache(::Type{IndicatorHennemannGassner},
1010
equations::AbstractEquations{2}, basis::LobattoLegendreBasis)
11-
alpha = Vector{real(basis)}()
11+
uEltype = real(basis)
12+
alpha = Vector{uEltype}()
1213
alpha_tmp = similar(alpha)
1314

14-
A = Array{real(basis), ndims(equations)}
15-
indicator_threaded = [A(undef, nnodes(basis), nnodes(basis))
16-
for _ in 1:Threads.maxthreadid()]
17-
modal_threaded = [A(undef, nnodes(basis), nnodes(basis))
18-
for _ in 1:Threads.maxthreadid()]
19-
modal_tmp1_threaded = [A(undef, nnodes(basis), nnodes(basis))
20-
for _ in 1:Threads.maxthreadid()]
15+
MA2d = MArray{Tuple{nnodes(basis), nnodes(basis)},
16+
uEltype, 2, nnodes(basis)^ndims(equations)}
17+
18+
indicator_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
19+
modal_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
20+
modal_tmp1_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
2121

2222
return (; alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded)
2323
end
@@ -127,13 +127,15 @@ function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}},
127127
end
128128

129129
# this method is used when the indicator is constructed as for shock-capturing volume integrals
130-
function create_cache(::Type{IndicatorLöhner}, equations::AbstractEquations{2},
131-
basis::LobattoLegendreBasis)
132-
alpha = Vector{real(basis)}()
130+
function create_cache(::Union{Type{IndicatorLöhner}, Type{IndicatorMax}},
131+
equations::AbstractEquations{2}, basis::LobattoLegendreBasis)
132+
uEltype = real(basis)
133+
alpha = Vector{uEltype}()
134+
135+
MA2d = MArray{Tuple{nnodes(basis), nnodes(basis)},
136+
uEltype, 2, nnodes(basis)^ndims(equations)}
133137

134-
A = Array{real(basis), ndims(equations)}
135-
indicator_threaded = [A(undef, nnodes(basis), nnodes(basis))
136-
for _ in 1:Threads.maxthreadid()]
138+
indicator_threaded = MA2d[MA2d(undef) for _ in 1:Threads.maxthreadid()]
137139

138140
return (; alpha, indicator_threaded)
139141
end
@@ -179,18 +181,6 @@ function (löhner::IndicatorLöhner)(u::AbstractArray{<:Any, 4},
179181
return alpha
180182
end
181183

182-
# this method is used when the indicator is constructed as for shock-capturing volume integrals
183-
function create_cache(::Type{IndicatorMax}, equations::AbstractEquations{2},
184-
basis::LobattoLegendreBasis)
185-
alpha = Vector{real(basis)}()
186-
187-
A = Array{real(basis), ndims(equations)}
188-
indicator_threaded = [A(undef, nnodes(basis), nnodes(basis))
189-
for _ in 1:Threads.maxthreadid()]
190-
191-
return (; alpha, indicator_threaded)
192-
end
193-
194184
function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 4},
195185
mesh, equations, dg::DGSEM, cache;
196186
kwargs...)

0 commit comments

Comments
 (0)