Skip to content

Commit 5bdb3c0

Browse files
DanielDoehringDanielDoehring
authored
Merge branch 'main' into split_Jac_SCT2
2 parents 75af937 + 6e51e7e commit 5bdb3c0

File tree

11 files changed

+180
-198
lines changed

11 files changed

+180
-198
lines changed

‎examples/p4est_2d_dgsem/elixir_navierstokes_blast_reflective.jl‎

Lines changed: 1 addition & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -94,10 +94,6 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
9494

9595
alive_callback = AliveCallback(analysis_interval = analysis_interval)
9696

97-
callbacks = CallbackSet(summary_callback,
98-
analysis_callback,
99-
alive_callback)
100-
10197
callbacks = CallbackSet(summary_callback,
10298
analysis_callback,
10399
alive_callback)
@@ -108,7 +104,5 @@ callbacks = CallbackSet(summary_callback,
108104
# https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Low-Storage-Methods
109105
ode_alg = CKLLSRK65_4M_4R()
110106

111-
abstol = 1e-6
112-
reltol = 1e-4
113-
sol = solve(ode, ode_alg; abstol = abstol, reltol = reltol,
107+
sol = solve(ode, ode_alg; abstol = 1e-6, reltol = 1e-4,
114108
ode_default_options()..., callback = callbacks);

‎examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -66,12 +66,15 @@ u_ls = A_matrix \ b
6666
# Iterative solve, works directly on the linear map, no explicit matrix construction required!
6767
using Krylov
6868

69+
atol = 1e-11
70+
rtol = 1e-10
71+
6972
# This solves the Laplace equation (i.e., steady-state diffusion/heat equation).
7073
# Note that we do not use CG since the operator `A_map` itself is only
7174
# symmetric and positive definite with respect to the inner product induced by
7275
# the DGSEM quadrature, not the standard Euclidean inner product. Standard CG
7376
# would require multiplying the result of `A_map * u` (and `b`) by the mass matrix.
74-
u_ls, stats = gmres(A_map, b)
77+
u_ls, stats = gmres(A_map, b, atol = atol, rtol = rtol)
7578

7679
###############################################################################
7780

‎src/auxiliary/containers.jl‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -350,6 +350,34 @@ function storage_type(C::Type{<:AbstractContainer})
350350
return storage_type(Adapt.unwrap_type(C))
351351
end
352352

353+
abstract type AbstractTreeElementContainer <: AbstractContainer end
354+
355+
# Return number of elements
356+
@inline nelements(elements::AbstractTreeElementContainer) = length(elements.cell_ids)
357+
# TODO: Taal performance, 1:nelements(elements) vs. Base.OneTo(nelements(elements))
358+
"""
359+
eachelement(elements::AbstractTreeElementContainer)
360+
361+
Return an iterator over the indices that specify the location in relevant data structures
362+
for the elements in `elements`.
363+
In particular, not the elements themselves are returned.
364+
"""
365+
@inline eachelement(elements::AbstractTreeElementContainer) = Base.OneTo(nelements(elements))
366+
@inline Base.real(elements::AbstractTreeElementContainer) = eltype(elements.node_coordinates)
367+
@inline nvariables(elements::AbstractTreeElementContainer) = size(elements.surface_flux_values,
368+
1)
369+
@inline nnodes(elements::AbstractTreeElementContainer) = size(elements.node_coordinates,
370+
2)
371+
@inline Base.eltype(elements::AbstractTreeElementContainer) = eltype(elements.surface_flux_values)
372+
373+
abstract type AbstractTreeBoundaryContainer <: AbstractContainer end
374+
375+
@inline nvariables(boundaries::AbstractTreeBoundaryContainer) = size(boundaries.u, 2)
376+
# Return number of boundaries
377+
@inline nboundaries(boundaries::AbstractTreeBoundaryContainer) = length(boundaries.orientations)
378+
# For 2D and 3D. 1D Hard-coded to 1
379+
@inline nnodes(boundaries::AbstractTreeBoundaryContainer) = size(boundaries.u, 3)
380+
353381
# backend handling
354382
"""
355383
trixi_backend(x)

‎src/auxiliary/precompile.jl‎

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -397,54 +397,54 @@ function _precompile_manual_()
397397
# 1D, serial
398398
@assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
399399
TreeMesh{1, Trixi.SerialTree{1}, RealT},
400-
Trixi.ElementContainer1D{RealT, uEltype}})
400+
Trixi.TreeElementContainer1D{RealT, uEltype}})
401401
@assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
402402
TreeMesh{1, Trixi.SerialTree{1}, RealT},
403-
Trixi.ElementContainer1D{RealT, uEltype}})
403+
Trixi.TreeElementContainer1D{RealT, uEltype}})
404404
@assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
405405
TreeMesh{1, Trixi.SerialTree{1}, RealT}, String})
406406

407407
# 2D, serial
408408
@assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
409409
TreeMesh{2, Trixi.SerialTree{2}, RealT},
410-
Trixi.ElementContainer2D{RealT, uEltype}})
410+
Trixi.TreeElementContainer2D{RealT, uEltype}})
411411
@assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
412412
TreeMesh{2, Trixi.SerialTree{2}, RealT},
413-
Trixi.ElementContainer2D{RealT, uEltype}})
413+
Trixi.TreeElementContainer2D{RealT, uEltype}})
414414
@assert Base.precompile(Tuple{typeof(Trixi.init_mortars), Array{Int, 1},
415415
TreeMesh{2, Trixi.SerialTree{2}, RealT},
416-
Trixi.ElementContainer2D{RealT, uEltype},
416+
Trixi.TreeElementContainer2D{RealT, uEltype},
417417
mortar_type})
418418
@assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
419419
TreeMesh{2, Trixi.SerialTree{2}, RealT}, String})
420420

421421
# 2D, parallel
422422
@assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
423423
TreeMesh{2, Trixi.ParallelTree{2}, RealT},
424-
Trixi.ElementContainer2D{RealT, uEltype}})
424+
Trixi.TreeElementContainer2D{RealT, uEltype}})
425425
@assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
426426
TreeMesh{2, Trixi.ParallelTree{2}, RealT},
427-
Trixi.ElementContainer2D{RealT, uEltype}})
427+
Trixi.TreeElementContainer2D{RealT, uEltype}})
428428
@assert Base.precompile(Tuple{typeof(Trixi.init_mortars), Array{Int, 1},
429429
TreeMesh{2, Trixi.ParallelTree{2}, RealT},
430-
Trixi.ElementContainer2D{RealT, uEltype},
430+
Trixi.TreeElementContainer2D{RealT, uEltype},
431431
mortar_type})
432432
@assert Base.precompile(Tuple{typeof(Trixi.init_mpi_interfaces), Array{Int, 1},
433433
TreeMesh{2, Trixi.ParallelTree{2}, RealT},
434-
Trixi.ElementContainer2D{RealT, uEltype}})
434+
Trixi.TreeElementContainer2D{RealT, uEltype}})
435435
@assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
436436
TreeMesh{2, Trixi.ParallelTree{2}, RealT}, String})
437437

438438
# 3D, serial
439439
@assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
440440
TreeMesh{3, Trixi.SerialTree{3}, RealT},
441-
Trixi.ElementContainer3D{RealT, uEltype}})
441+
Trixi.TreeElementContainer3D{RealT, uEltype}})
442442
@assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
443443
TreeMesh{3, Trixi.SerialTree{3}, RealT},
444-
Trixi.ElementContainer3D{RealT, uEltype}})
444+
Trixi.TreeElementContainer3D{RealT, uEltype}})
445445
@assert Base.precompile(Tuple{typeof(Trixi.init_mortars), Array{Int, 1},
446446
TreeMesh{3, Trixi.SerialTree{3}, RealT},
447-
Trixi.ElementContainer3D{RealT, uEltype},
447+
Trixi.TreeElementContainer3D{RealT, uEltype},
448448
mortar_type})
449449
@assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
450450
TreeMesh{3, Trixi.SerialTree{3}, RealT}, String})

‎src/solvers/dgsem_structured/containers.jl‎

Lines changed: 15 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -5,8 +5,8 @@
55
@muladd begin
66
#! format: noindent
77

8-
struct ElementContainer{NDIMS, RealT <: Real, uEltype <: Real,
9-
NDIMSP1, NDIMSP2, NDIMSP3}
8+
struct StructuredElementContainer{NDIMS, RealT <: Real, uEltype <: Real,
9+
NDIMSP1, NDIMSP2, NDIMSP3}
1010
# Physical coordinates at each node
1111
node_coordinates::Array{RealT, NDIMSP2} # [orientation, node_i, node_j, node_k, element]
1212
# ID of neighbor element in negative direction in orientation
@@ -49,22 +49,25 @@ function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT},
4949
NDIMS - 1)..., NDIMS * 2,
5050
nelements)
5151

52-
elements = ElementContainer{NDIMS, RealT, uEltype, NDIMS + 1, NDIMS + 2, NDIMS + 3}(node_coordinates,
53-
left_neighbors,
54-
jacobian_matrix,
55-
contravariant_vectors,
56-
inverse_jacobian,
57-
surface_flux_values)
52+
elements = StructuredElementContainer{NDIMS, RealT, uEltype,
53+
NDIMS + 1, NDIMS + 2, NDIMS + 3}(node_coordinates,
54+
left_neighbors,
55+
jacobian_matrix,
56+
contravariant_vectors,
57+
inverse_jacobian,
58+
surface_flux_values)
5859

5960
init_elements!(elements, mesh, basis)
6061
return elements
6162
end
6263

63-
@inline nelements(elements::ElementContainer) = size(elements.left_neighbors, 2)
64-
@inline Base.ndims(::ElementContainer{NDIMS}) where {NDIMS} = NDIMS
64+
@inline nelements(elements::StructuredElementContainer) = size(elements.left_neighbors,
65+
2)
6566

66-
function Base.eltype(::ElementContainer{NDIMS, RealT, uEltype}) where {NDIMS, RealT,
67-
uEltype}
67+
function Base.eltype(::StructuredElementContainer{NDIMS, RealT, uEltype}) where {NDIMS,
68+
RealT,
69+
uEltype
70+
}
6871
return uEltype
6972
end
7073

‎src/solvers/dgsem_tree/containers_1d.jl‎

Lines changed: 29 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,8 @@
66
#! format: noindent
77

88
# Container data structure (structure-of-arrays style) for DG elements
9-
mutable struct ElementContainer1D{RealT <: Real, uEltype <: Real} <: AbstractContainer
9+
mutable struct TreeElementContainer1D{RealT <: Real, uEltype <: Real} <:
10+
AbstractTreeElementContainer
1011
inverse_jacobian::Vector{RealT} # [elements]
1112
node_coordinates::Array{RealT, 3} # [orientation, i, elements]
1213
surface_flux_values::Array{uEltype, 3} # [variables, direction, elements]
@@ -16,16 +17,12 @@ mutable struct ElementContainer1D{RealT <: Real, uEltype <: Real} <: AbstractCon
1617
_surface_flux_values::Vector{uEltype}
1718
end
1819

19-
nvariables(elements::ElementContainer1D) = size(elements.surface_flux_values, 1)
20-
nnodes(elements::ElementContainer1D) = size(elements.node_coordinates, 2)
21-
Base.eltype(elements::ElementContainer1D) = eltype(elements.surface_flux_values)
22-
2320
# Only one-dimensional `Array`s are `resize!`able in Julia.
2421
# Hence, we use `Vector`s as internal storage and `resize!`
2522
# them whenever needed. Then, we reuse the same memory by
2623
# `unsafe_wrap`ping multi-dimensional `Array`s around the
2724
# internal storage.
28-
function Base.resize!(elements::ElementContainer1D, capacity)
25+
function Base.resize!(elements::TreeElementContainer1D, capacity)
2926
n_nodes = nnodes(elements)
3027
n_variables = nvariables(elements)
3128
@unpack _node_coordinates, _surface_flux_values,
@@ -46,9 +43,9 @@ function Base.resize!(elements::ElementContainer1D, capacity)
4643
return nothing
4744
end
4845

49-
function ElementContainer1D{RealT, uEltype}(capacity::Integer, n_variables,
50-
n_nodes) where {RealT <: Real,
51-
uEltype <: Real}
46+
function TreeElementContainer1D{RealT, uEltype}(capacity::Integer, n_variables,
47+
n_nodes) where {RealT <: Real,
48+
uEltype <: Real}
5249
nan_RealT = convert(RealT, NaN)
5350
nan_uEltype = convert(uEltype, NaN)
5451

@@ -65,33 +62,21 @@ function ElementContainer1D{RealT, uEltype}(capacity::Integer, n_variables,
6562

6663
cell_ids = fill(typemin(Int), capacity)
6764

68-
return ElementContainer1D{RealT, uEltype}(inverse_jacobian, node_coordinates,
69-
surface_flux_values, cell_ids,
70-
_node_coordinates, _surface_flux_values)
65+
return TreeElementContainer1D{RealT, uEltype}(inverse_jacobian, node_coordinates,
66+
surface_flux_values, cell_ids,
67+
_node_coordinates,
68+
_surface_flux_values)
7169
end
7270

73-
# Return number of elements
74-
@inline nelements(elements::ElementContainer1D) = length(elements.cell_ids)
75-
# TODO: Taal performance, 1:nelements(elements) vs. Base.OneTo(nelements(elements))
76-
"""
77-
eachelement(elements::ElementContainer1D)
78-
79-
Return an iterator over the indices that specify the location in relevant data structures
80-
for the elements in `elements`.
81-
In particular, not the elements themselves are returned.
82-
"""
83-
@inline eachelement(elements::ElementContainer1D) = Base.OneTo(nelements(elements))
84-
@inline Base.real(elements::ElementContainer1D) = eltype(elements.node_coordinates)
85-
8671
# Create element container and initialize element data
8772
function init_elements(cell_ids, mesh::TreeMesh1D,
8873
equations::AbstractEquations{1},
8974
basis, ::Type{RealT},
9075
::Type{uEltype}) where {RealT <: Real, uEltype <: Real}
9176
# Initialize container
9277
n_elements = length(cell_ids)
93-
elements = ElementContainer1D{RealT, uEltype}(n_elements, nvariables(equations),
94-
nnodes(basis))
78+
elements = TreeElementContainer1D{RealT, uEltype}(n_elements, nvariables(equations),
79+
nnodes(basis))
9580

9681
init_elements!(elements, cell_ids, mesh, basis)
9782
return elements
@@ -192,7 +177,7 @@ end
192177

193178
# Create interface container and initialize interface data in `elements`.
194179
function init_interfaces(cell_ids, mesh::TreeMesh1D,
195-
elements::ElementContainer1D)
180+
elements::TreeElementContainer1D)
196181
# Initialize container
197182
n_interfaces = count_required_interfaces(mesh, cell_ids)
198183
interfaces = InterfaceContainer1D{eltype(elements)}(n_interfaces,
@@ -284,7 +269,8 @@ function init_interfaces!(interfaces, elements, mesh::TreeMesh1D)
284269
end
285270

286271
# Container data structure (structure-of-arrays style) for DG boundaries
287-
mutable struct BoundaryContainer1D{RealT <: Real, uEltype <: Real} <: AbstractContainer
272+
mutable struct TreeBoundaryContainer1D{RealT <: Real, uEltype <: Real} <:
273+
AbstractTreeBoundaryContainer
288274
u::Array{uEltype, 3} # [leftright, variables, boundaries]
289275
neighbor_ids::Vector{Int} # [boundaries]
290276
orientations::Vector{Int} # [boundaries]
@@ -296,11 +282,11 @@ mutable struct BoundaryContainer1D{RealT <: Real, uEltype <: Real} <: AbstractCo
296282
_node_coordinates::Vector{RealT}
297283
end
298284

299-
nvariables(boundaries::BoundaryContainer1D) = size(boundaries.u, 2)
300-
Base.eltype(boundaries::BoundaryContainer1D) = eltype(boundaries.u)
285+
# 1D: Only one boundary node
286+
nnodes(boundaries::TreeBoundaryContainer1D) = 1
301287

302288
# See explanation of Base.resize! for the element container
303-
function Base.resize!(boundaries::BoundaryContainer1D, capacity)
289+
function Base.resize!(boundaries::TreeBoundaryContainer1D, capacity)
304290
n_variables = nvariables(boundaries)
305291
@unpack _u, _node_coordinates,
306292
neighbor_ids, orientations, neighbor_sides = boundaries
@@ -322,9 +308,9 @@ function Base.resize!(boundaries::BoundaryContainer1D, capacity)
322308
return nothing
323309
end
324310

325-
function BoundaryContainer1D{RealT, uEltype}(capacity::Integer, n_variables,
326-
n_nodes) where {RealT <: Real,
327-
uEltype <: Real}
311+
function TreeBoundaryContainer1D{RealT, uEltype}(capacity::Integer,
312+
n_variables) where {RealT <: Real,
313+
uEltype <: Real}
328314
nan_RealT = convert(RealT, NaN)
329315
nan_uEltype = convert(uEltype, NaN)
330316

@@ -345,24 +331,20 @@ function BoundaryContainer1D{RealT, uEltype}(capacity::Integer, n_variables,
345331

346332
n_boundaries_per_direction = SVector(0, 0)
347333

348-
return BoundaryContainer1D{RealT, uEltype}(u, neighbor_ids, orientations,
349-
neighbor_sides,
350-
node_coordinates,
351-
n_boundaries_per_direction,
352-
_u, _node_coordinates)
334+
return TreeBoundaryContainer1D{RealT, uEltype}(u, neighbor_ids, orientations,
335+
neighbor_sides,
336+
node_coordinates,
337+
n_boundaries_per_direction,
338+
_u, _node_coordinates)
353339
end
354340

355-
# Return number of boundaries
356-
nboundaries(boundaries::BoundaryContainer1D) = length(boundaries.orientations)
357-
358341
# Create boundaries container and initialize boundary data in `elements`.
359342
function init_boundaries(cell_ids, mesh::TreeMesh1D,
360-
elements::ElementContainer1D)
343+
elements::TreeElementContainer1D)
361344
# Initialize container
362345
n_boundaries = count_required_boundaries(mesh, cell_ids)
363-
boundaries = BoundaryContainer1D{real(elements), eltype(elements)}(n_boundaries,
364-
nvariables(elements),
365-
nnodes(elements))
346+
boundaries = TreeBoundaryContainer1D{real(elements), eltype(elements)}(n_boundaries,
347+
nvariables(elements))
366348

367349
# Connect elements with boundaries
368350
init_boundaries!(boundaries, elements, mesh)

0 commit comments

Comments
 (0)