26328518171

Committed 23 May 2026 08:49AM UTC coverage: 94.161% (-2.7%) from 96.896%

Build # 26328518171

Build Type

push

github

Committed by

web-flow

Commit Message

hotfix saving restart files using GPUs (#3032)

Coverage Stats

3 of 3 new or added lines in 1 file covered. (100.0%)

1346 existing lines in 63 files now uncovered.

46350 of 49224 relevant lines covered (94.16%)

39984412.7 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

19.4

/src/postprocessing/spectral_analysis_2d.jl

# By default, Julia/LLVM does not use fused multiply-add operations (FMAs)
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details
@muladd begin
#! format: noindent

"""
    compute_kinetic_energy_spectrum(v1, v2)

Compute an isotropic 1D kinetic energy spectrum from two 2D Cartesian velocity
components `v1` and `v2`. For compressible Euler kinetic energy spectra,
pass density-weighted components `sqrt(rho) * v1` and `sqrt(rho) * v2`.
The modal energy is normalized by `1 / N^2'.
"""
function compute_kinetic_energy_spectrum(v1::AbstractArray{<:Any, 2},
                                         v2::AbstractArray{<:Any, 2})

    # Compute the energy modes using FFTW
    energy_modes = 0.5f0 .* (abs2.(fft(v1)) .+ abs2.(fft(v2)))
    energy_modes ./= length(energy_modes)^2

    return radial_energy_spectrum(energy_modes)
end

"""
    compute_kinetic_energy_spectrum(u, mesh::TreeMesh{2}, equations, solver::DGSEM,
                                    cache)

Compute the energy spectrum for a non-AMR 2D `TreeMesh`/`DGSEM` solution by first
interpolating from LGL nodes to a uniform Cartesian grid.
"""
function compute_kinetic_energy_spectrum(u, mesh::TreeMesh{2},
                                         equations::AbstractCompressibleEulerEquations,
                                         solver::DGSEM, cache)
    # Interpolates conservative polynomials to a uniform Cartesian grid then converts to primitives at each uniform node
    u_uniform = interpolate_lgl_to_uniform_cartesian(u, mesh, equations, solver, cache)
    grid_size = size(u_uniform)[2:end] # the first dimension is the equation index so it is not needed to count the spatial indices
    rho = Array{eltype(u)}(undef, grid_size)
    v1 = Array{eltype(u)}(undef, grid_size)
    v2 = Array{eltype(u)}(undef, grid_size)
    for idx in CartesianIndices(grid_size)
        u_node = get_node_vars(u_uniform, equations, solver, Tuple(idx)...)
        prim = cons2prim(u_node, equations)
        rho[idx] = prim[1]
        v1[idx] = prim[2]
        v2[idx] = prim[3]
    end
    # Convert primitive velocity components to density weighted form before FFT
    density_weighted_velocity_1 = sqrt.(rho) .* v1
    density_weighted_velocity_2 = sqrt.(rho) .* v2

    return compute_kinetic_energy_spectrum(density_weighted_velocity_1,
                                           density_weighted_velocity_2)
end

function interpolate_lgl_to_uniform_cartesian(u, mesh::TreeMesh{2},
                                              equations::AbstractCompressibleEulerEquations,
                                              solver::DGSEM, cache)
    # Restrict to straightforward non AMR setups
    leaf_cell_ids = leaf_cells(mesh.tree)
    levels = mesh.tree.levels[leaf_cell_ids]
    if !all(==(first(levels)), levels)
        throw(ArgumentError("Non-uniform meshes are not supported yet"))
    end

    level = first(levels)
    cells_per_dimension = 2^level

    # Uses one uniform interpolation node per DGSEM solution node in each coordinate
    # direction. A degree p element has p + 1 LGL nodes per direction, so this
    # is the minimum Cartesian sampling matching the element-wise geometry
    # but this can be increased to finer sampling sizes if needed
    n_uniform_nodes = polydeg(solver) + 1
    grid_points_per_dimension = n_uniform_nodes * cells_per_dimension

    n_vars = nvariables(equations)
    u_uniform = Array{eltype(u)}(undef, n_vars, grid_points_per_dimension,
                                 grid_points_per_dimension)

    # Interpolate from LGL nodes to cell-centered equidistant nodes in each element
    dx_reference = 2 / n_uniform_nodes
    nodes_out = collect(range(-1 + dx_reference / 2, 1 - dx_reference / 2,
                              length = n_uniform_nodes))
    vandermonde = polynomial_interpolation_matrix(get_nodes(solver.basis), nodes_out)

    # Element center coordinates determine where each interpolated block sits in the global grid
    center = reshape(mesh.tree.center_level_0, :, 1)
    normalized_coordinates = (mesh.tree.coordinates[:, leaf_cell_ids] .- center) ./
                             (mesh.tree.length_level_0 / 2)
    dx_global = 2 / (n_uniform_nodes * cells_per_dimension)

    # Example for 1st order (2 nodes per direction per element) for a 2x2 tree of cells which produces a 4x4 global FFT grid.
    #
    # Reference element — local indices (i along ξ, j along η):
    #
    #      j=2  (1,2)--------(2,2)
    #      ↑     |             |
    #      |     |             |
    #      |     |             |
    #      |     |             |
    #      j    (1,1)--------(2,1)
    #             ---------> i (ξ)
    #
    #   Each of the 4 tree cells pastes a 2×2 block into the 4×4 array global grid, 'u_uniform'.
    #   Labels (column, row) are indices into that array: column ↔ x, row ↔ y on the domain.
    #   The assignment `u_uniform[:, r1, r2] .= interpolated` fills one such rectangle.
    #   The bottom-right element of the original 2×2 block is labeled as 'E' in the diagram below and fills
    #   rows 3 to 4 and columns 3 to 4 of the global grid.
    #
    #         column  1         2         3         4
    #                 +---------+---------+---------+---------+
    #         row 1   | (1,1)   | (2,1)   | (3,1)   | (4,1)   |
    #                 +---------+---------+---------+---------+
    #         row 2   | (1,2)   | (2,2)   | (3,2)   | (4,2)   |
    #                 +---------+---------+---------+---------+   ↑
    #         row 3   | (1,3)   | (2,3)   | (3,3) E | (4,3) E |   |
    #                 +---------+---------+---------+---------+   |  r2 = 3:4
    #         row 4   | (1,4)   | (2,4)   | (3,4) E | (4,4) E |   |
    #                 +---------+---------+---------+---------+   ↓
    #                                     |<-----r1 = 3:4---->|
    #

    for element in eachelement(solver, cache)
        # Gather conservative nodal values on the reference LGL tensor grid for the element
        u_sample = get_node_vars(u, equations, solver, 1, 1, element)
        element_size = (n_vars, nnodes(solver), nnodes(solver))
        element_values = Array{eltype(u_sample)}(undef,
                                                 element_size)
        for j in eachnode(solver), i in eachnode(solver)
            u_node = get_node_vars(u, equations, solver, i, j, element)
            for variable in 1:n_vars
                element_values[variable, i, j] = u_node[variable]
            end
        end
        interpolated = multiply_dimensionwise(vandermonde, element_values)

        # Each element is placed on the uniform grid assuming reference directions align with
        # physical axes (ξ→x, η→y) and nodal values use the DGSEM tensor order along (ξ, η).
        # first_index[dim] refers to the first index of the nodes of this element within the array of global tensor nodes
        first_index = Vector{Int}(undef, 2)
        for dim in 1:2
            lower_left = normalized_coordinates[dim, element] -
                         (n_uniform_nodes - 1) / 2 * dx_global
            first_index[dim] = round(Int,
                                     (lower_left - (-1 + dx_global / 2)) /
                                     dx_global) + 1
        end

        # with `u_uniform` as the full uniform Cartesian grid over the domain, `r1` and `r2` are the
        # global grid indices corresponding to `u_uniform` that this specific element's interpolated block fits within
        # Essentially, `r1` and `r2` are the positions of the current element within `u_uniform`
        r1 = first_index[1]:(first_index[1] + n_uniform_nodes - 1)
        r2 = first_index[2]:(first_index[2] + n_uniform_nodes - 1)
        u_uniform[:, r1, r2] .= interpolated[:, :, :]
    end
    return u_uniform
end

"""
    compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{2}, equations,
                                    dg::DGMultiSBP, cache)

Compute the energy spectrum for a 2D `DGMulti` finite-difference SBP solution whose
nodes already form a uniform Cartesian grid.
"""
function compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{2},
                                         equations::AbstractCompressibleEulerEquations,
                                         dg::DGMultiSBP, cache)
    # Unpacks the primitive variables from the conservative state for FDSBP DGMulti solutions
    u_values = parent(u)
    n_points = length(u_values)
    n = round(Int, sqrt(n_points))
    q = cons2prim.(u_values, Ref(equations)) # q is the vector that contains the primitive variables for density and velocity converted from the conservative variables
    rho = reshape(getindex.(q, 1), n, n)
    density_weighted_velocity_1 = sqrt.(rho) .* reshape(getindex.(q, 2), n, n)
    density_weighted_velocity_2 = sqrt.(rho) .* reshape(getindex.(q, 3), n, n)

    return compute_kinetic_energy_spectrum(density_weighted_velocity_1,
                                           density_weighted_velocity_2)
end
end # @muladd

1	# By default, Julia/LLVM does not use fused multiply-add operations (FMAs)
2	# Since these FMAs can increase the performance of many numerical algorithms,
3	# we need to opt-in explicitly
4	# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details
5	@muladd begin
6	#! format: noindent
7
8	"""
9	compute_kinetic_energy_spectrum(v1, v2)
10
11	Compute an isotropic 1D kinetic energy spectrum from two 2D Cartesian velocity
12	components `v1` and `v2`. For compressible Euler kinetic energy spectra,
13	pass density-weighted components `sqrt(rho) * v1` and `sqrt(rho) * v2`.
14	The modal energy is normalized by `1 / N^2'.
15	"""
16	function compute_kinetic_energy_spectrum(v1::AbstractArray{<:Any, 2},	3✔
17	v2::AbstractArray{<:Any, 2})
18
19	# Compute the energy modes using FFTW
20	energy_modes = 0.5f0 .* (abs2.(fft(v1)) .+ abs2.(fft(v2)))	6✔
21	energy_modes ./= length(energy_modes)^2	6✔
22
23	return radial_energy_spectrum(energy_modes)	3✔
24	end
25
26	"""
27	compute_kinetic_energy_spectrum(u, mesh::TreeMesh{2}, equations, solver::DGSEM,
28	cache)
29
30	Compute the energy spectrum for a non-AMR 2D `TreeMesh`/`DGSEM` solution by first
31	interpolating from LGL nodes to a uniform Cartesian grid.
32	"""
UNCOV 33	function compute_kinetic_energy_spectrum(u, mesh::TreeMesh{2},	×
34	equations::AbstractCompressibleEulerEquations,
35	solver::DGSEM, cache)
36	# Interpolates conservative polynomials to a uniform Cartesian grid then converts to primitives at each uniform node
UNCOV 37	u_uniform = interpolate_lgl_to_uniform_cartesian(u, mesh, equations, solver, cache)	×
UNCOV 38	grid_size = size(u_uniform)[2:end] # the first dimension is the equation index so it is not needed to count the spatial indices	×
UNCOV 39	rho = Array{eltype(u)}(undef, grid_size)	×
UNCOV 40	v1 = Array{eltype(u)}(undef, grid_size)	×
UNCOV 41	v2 = Array{eltype(u)}(undef, grid_size)	×
UNCOV 42	for idx in CartesianIndices(grid_size)	×
UNCOV 43	u_node = get_node_vars(u_uniform, equations, solver, Tuple(idx)...)	×
UNCOV 44	prim = cons2prim(u_node, equations)	×
UNCOV 45	rho[idx] = prim[1]	×
UNCOV 46	v1[idx] = prim[2]	×
UNCOV 47	v2[idx] = prim[3]	×
UNCOV 48	end	×
49	# Convert primitive velocity components to density weighted form before FFT
UNCOV 50	density_weighted_velocity_1 = sqrt.(rho) .* v1	×
UNCOV 51	density_weighted_velocity_2 = sqrt.(rho) .* v2	×
52
UNCOV 53	return compute_kinetic_energy_spectrum(density_weighted_velocity_1,	×
54	density_weighted_velocity_2)
55	end
56
UNCOV 57	function interpolate_lgl_to_uniform_cartesian(u, mesh::TreeMesh{2},	×
58	equations::AbstractCompressibleEulerEquations,
59	solver::DGSEM, cache)
60	# Restrict to straightforward non AMR setups
UNCOV 61	leaf_cell_ids = leaf_cells(mesh.tree)	×
UNCOV 62	levels = mesh.tree.levels[leaf_cell_ids]	×
UNCOV 63	if !all(==(first(levels)), levels)	×
64	throw(ArgumentError("Non-uniform meshes are not supported yet"))	×
65	end
66
UNCOV 67	level = first(levels)	×
UNCOV 68	cells_per_dimension = 2^level	×
69
70	# Uses one uniform interpolation node per DGSEM solution node in each coordinate
71	# direction. A degree p element has p + 1 LGL nodes per direction, so this
72	# is the minimum Cartesian sampling matching the element-wise geometry
73	# but this can be increased to finer sampling sizes if needed
UNCOV 74	n_uniform_nodes = polydeg(solver) + 1	×
UNCOV 75	grid_points_per_dimension = n_uniform_nodes * cells_per_dimension	×
76
UNCOV 77	n_vars = nvariables(equations)	×
UNCOV 78	u_uniform = Array{eltype(u)}(undef, n_vars, grid_points_per_dimension,	×
79	grid_points_per_dimension)
80
81	# Interpolate from LGL nodes to cell-centered equidistant nodes in each element
UNCOV 82	dx_reference = 2 / n_uniform_nodes	×
UNCOV 83	nodes_out = collect(range(-1 + dx_reference / 2, 1 - dx_reference / 2,	×
84	length = n_uniform_nodes))
UNCOV 85	vandermonde = polynomial_interpolation_matrix(get_nodes(solver.basis), nodes_out)	×
86
87	# Element center coordinates determine where each interpolated block sits in the global grid
UNCOV 88	center = reshape(mesh.tree.center_level_0, :, 1)	×
UNCOV 89	normalized_coordinates = (mesh.tree.coordinates[:, leaf_cell_ids] .- center) ./	×
90	(mesh.tree.length_level_0 / 2)
UNCOV 91	dx_global = 2 / (n_uniform_nodes * cells_per_dimension)	×
92
93	# Example for 1st order (2 nodes per direction per element) for a 2x2 tree of cells which produces a 4x4 global FFT grid.
94	#
95	# Reference element — local indices (i along ξ, j along η):
96	#
97	# j=2 (1,2)--------(2,2)
98	# ↑ \| \|
99	# \| \| \|
100	# \| \| \|
101	# \| \| \|
102	# j (1,1)--------(2,1)
103	# ---------> i (ξ)
104	#
105	# Each of the 4 tree cells pastes a 2×2 block into the 4×4 array global grid, 'u_uniform'.
106	# Labels (column, row) are indices into that array: column ↔ x, row ↔ y on the domain.
107	# The assignment `u_uniform[:, r1, r2] .= interpolated` fills one such rectangle.
108	# The bottom-right element of the original 2×2 block is labeled as 'E' in the diagram below and fills
109	# rows 3 to 4 and columns 3 to 4 of the global grid.
110	#
111	# column 1 2 3 4
112	# +---------+---------+---------+---------+
113	# row 1 \| (1,1) \| (2,1) \| (3,1) \| (4,1) \|
114	# +---------+---------+---------+---------+
115	# row 2 \| (1,2) \| (2,2) \| (3,2) \| (4,2) \|
116	# +---------+---------+---------+---------+ ↑
117	# row 3 \| (1,3) \| (2,3) \| (3,3) E \| (4,3) E \| \|
118	# +---------+---------+---------+---------+ \| r2 = 3:4
119	# row 4 \| (1,4) \| (2,4) \| (3,4) E \| (4,4) E \| \|
120	# +---------+---------+---------+---------+ ↓
121	# \|<-----r1 = 3:4---->\|
122	#
123
UNCOV 124	for element in eachelement(solver, cache)	×
125	# Gather conservative nodal values on the reference LGL tensor grid for the element
UNCOV 126	u_sample = get_node_vars(u, equations, solver, 1, 1, element)	×
UNCOV 127	element_size = (n_vars, nnodes(solver), nnodes(solver))	×
UNCOV 128	element_values = Array{eltype(u_sample)}(undef,	×
129	element_size)
UNCOV 130	for j in eachnode(solver), i in eachnode(solver)	×
UNCOV 131	u_node = get_node_vars(u, equations, solver, i, j, element)	×
UNCOV 132	for variable in 1:n_vars	×
UNCOV 133	element_values[variable, i, j] = u_node[variable]	×
UNCOV 134	end	×
UNCOV 135	end	×
UNCOV 136	interpolated = multiply_dimensionwise(vandermonde, element_values)	×
137
138	# Each element is placed on the uniform grid assuming reference directions align with
139	# physical axes (ξ→x, η→y) and nodal values use the DGSEM tensor order along (ξ, η).
140	# first_index[dim] refers to the first index of the nodes of this element within the array of global tensor nodes
UNCOV 141	first_index = Vector{Int}(undef, 2)	×
UNCOV 142	for dim in 1:2	×
UNCOV 143	lower_left = normalized_coordinates[dim, element] -	×
144	(n_uniform_nodes - 1) / 2 * dx_global
UNCOV 145	first_index[dim] = round(Int,	×
146	(lower_left - (-1 + dx_global / 2)) /
147	dx_global) + 1
UNCOV 148	end	×
149
150	# with `u_uniform` as the full uniform Cartesian grid over the domain, `r1` and `r2` are the
151	# global grid indices corresponding to `u_uniform` that this specific element's interpolated block fits within
152	# Essentially, `r1` and `r2` are the positions of the current element within `u_uniform`
UNCOV 153	r1 = first_index[1]:(first_index[1] + n_uniform_nodes - 1)	×
UNCOV 154	r2 = first_index[2]:(first_index[2] + n_uniform_nodes - 1)	×
UNCOV 155	u_uniform[:, r1, r2] .= interpolated[:, :, :]	×
UNCOV 156	end	×
UNCOV 157	return u_uniform	×
158	end
159
160	"""
161	compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{2}, equations,
162	dg::DGMultiSBP, cache)
163
164	Compute the energy spectrum for a 2D `DGMulti` finite-difference SBP solution whose
165	nodes already form a uniform Cartesian grid.
166	"""
167	function compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{2},	1✔
168	equations::AbstractCompressibleEulerEquations,
169	dg::DGMultiSBP, cache)
170	# Unpacks the primitive variables from the conservative state for FDSBP DGMulti solutions
171	u_values = parent(u)	1✔
172	n_points = length(u_values)	1✔
173	n = round(Int, sqrt(n_points))	1✔
174	q = cons2prim.(u_values, Ref(equations)) # q is the vector that contains the primitive variables for density and velocity converted from the conservative variables	2✔
175	rho = reshape(getindex.(q, 1), n, n)	1✔
176	density_weighted_velocity_1 = sqrt.(rho) .* reshape(getindex.(q, 2), n, n)	2✔
177	density_weighted_velocity_2 = sqrt.(rho) .* reshape(getindex.(q, 3), n, n)	2✔
178
179	return compute_kinetic_energy_spectrum(density_weighted_velocity_1,	1✔
180	density_weighted_velocity_2)
181	end
182	end # @muladd

trixi-framework / Trixi.jl / 26328518171

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous