30296365372

Committed 27 Jul 2026 07:00PM UTC coverage: 85.072% (-5.9%) from 90.934%

Build # 30296365372

Build Type

Pull #3157

github

Committed by

web-flow

Commit Message

Merge 9443bad00 into b2ff19a99

Pull Request Pull Request #3157: fix `Trixi.Threaded`

Coverage Stats

43659 of 51320 relevant lines covered (85.07%)

32332827.23 hits per line

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

84.72

/src/postprocessing/spectral_analysis_3d.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, v3)

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

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

    return radial_energy_spectrum(energy_modes)
end

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

Compute the energy spectrum for a non-AMR 3D `TreeMesh`/`DGSEM` solution by first
interpolating from LGL nodes to a uniform Cartesian grid.
"""
function compute_kinetic_energy_spectrum(u, mesh::TreeMesh{3},
                                         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)
    v3 = 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]
        v3[idx] = prim[4]
    end
    # Converts primitive velocity components to density weighted form before FFT
    density_weighted_velocity_1 = sqrt.(rho) .* v1
    density_weighted_velocity_2 = sqrt.(rho) .* v2
    density_weighted_velocity_3 = sqrt.(rho) .* v3

    return compute_kinetic_energy_spectrum(density_weighted_velocity_1,
                                           density_weighted_velocity_2,
                                           density_weighted_velocity_3)
end

function interpolate_lgl_to_uniform_cartesian(u, mesh::TreeMesh{3},
                                              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,
                                 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)

    for element in eachelement(solver, cache)
        # Gather conservative nodal values on the reference LGL grid for the element
        u_sample = get_node_vars(u, equations, solver, 1, 1, 1, element)
        element_size = (n_vars, nnodes(solver), nnodes(solver),
                        nnodes(solver))
        element_values = Array{eltype(u_sample)}(undef, element_size)
        for k in eachnode(solver), j in eachnode(solver), i in eachnode(solver)
            u_node = get_node_vars(u, equations, solver, i, j, k, element)
            for variable in 1:n_vars
                element_values[variable, i, j, k] = 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, ζ→z) and nodal values use the standard DGSEM tensor order along (ξ, η, ζ)
        first_index = Vector{Int}(undef, 3)
        for dim in 1:3
            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

        # Writes the interpolated block onto the global grid for the larger output
        # `r1`, `r2`, and `r3` are the global indices corresponding to `u_uniform` that this specific element's interpolated block fits within
        # Essentially, `r1`, `r2`, and `r3` are the positions of the current element within `u_uniform`
        # first_index[dim] refers to the first index of the nodes of this element within the array of global tensor nodes
        # See the sketch in `spectral_analysis_2d.jl` for a sketch in 2d of how the local to global assembly works
        r1 = first_index[1]:(first_index[1] + n_uniform_nodes - 1)
        r2 = first_index[2]:(first_index[2] + n_uniform_nodes - 1)
        r3 = first_index[3]:(first_index[3] + n_uniform_nodes - 1)
        u_uniform[:, r1, r2, r3] .= interpolated[:, :, :, :]
    end
    return u_uniform
end

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

Compute the energy spectrum for a 3D `DGMulti` finite-difference SBP solution whose
nodes already form a uniform Cartesian grid.
"""
function compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{3},
                                         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, n_points^(1 / 3))
    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, n)
    density_weighted_velocity_1 = sqrt.(rho) .* reshape(getindex.(q, 2), n, n, n)
    density_weighted_velocity_2 = sqrt.(rho) .* reshape(getindex.(q, 3), n, n, n)
    density_weighted_velocity_3 = sqrt.(rho) .* reshape(getindex.(q, 4), n, n, n)

    return compute_kinetic_energy_spectrum(density_weighted_velocity_1,
                                           density_weighted_velocity_2,
                                           density_weighted_velocity_3)
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, v3)
10
11	Compute an isotropic 1D kinetic energy spectrum from three 3D Cartesian velocity
12	components `v1`, `v2`, `v3`. For compressible Euler kinetic energy spectra,
13	pass density-weighted components `sqrt(rho) * v1`, `sqrt(rho) * v2`, and `sqrt(rho) * v3`.
14	The modal energy is normalized by `1 / N^3`.
15	"""
16	function compute_kinetic_energy_spectrum(v1::AbstractArray{<:Any, 3},	3✔
17	v2::AbstractArray{<:Any, 3},
18	v3::AbstractArray{<:Any, 3})
19
20	# Compute the energy modes using FFTW
21	energy_modes = 0.5f0 .* (abs2.(fft(v1)) .+ abs2.(fft(v2)) .+ abs2.(fft(v3)))	6✔
22	energy_modes ./= length(energy_modes)^2	6✔
23
24	return radial_energy_spectrum(energy_modes)	3✔
25	end
26
27	"""
28	compute_kinetic_energy_spectrum(u, mesh::TreeMesh{3}, equations, solver::DGSEM,
29	cache)
30
31	Compute the energy spectrum for a non-AMR 3D `TreeMesh`/`DGSEM` solution by first
32	interpolating from LGL nodes to a uniform Cartesian grid.
33	"""
34	function compute_kinetic_energy_spectrum(u, mesh::TreeMesh{3},	1✔
35	equations::AbstractCompressibleEulerEquations,
36	solver::DGSEM, cache)
37	# Interpolates conservative polynomials to a uniform Cartesian grid then converts to primitives at each uniform node
38	u_uniform = interpolate_lgl_to_uniform_cartesian(u, mesh, equations, solver, cache)	1✔
39	grid_size = size(u_uniform)[2:end] # the first dimension is the equation index so it is not needed to count the spatial indices	1✔
40	rho = Array{eltype(u)}(undef, grid_size)	1✔
41	v1 = Array{eltype(u)}(undef, grid_size)	1✔
42	v2 = Array{eltype(u)}(undef, grid_size)	1✔
43	v3 = Array{eltype(u)}(undef, grid_size)	1✔
44	for idx in CartesianIndices(grid_size)	1✔
45	u_node = get_node_vars(u_uniform, equations, solver, Tuple(idx)...)	32,768✔
46	prim = cons2prim(u_node, equations)	32,768✔
47	rho[idx] = prim[1]	32,768✔
48	v1[idx] = prim[2]	32,768✔
49	v2[idx] = prim[3]	32,768✔
50	v3[idx] = prim[4]	32,768✔
51	end	32,769✔
52	# Converts primitive velocity components to density weighted form before FFT
53	density_weighted_velocity_1 = sqrt.(rho) .* v1	2✔
54	density_weighted_velocity_2 = sqrt.(rho) .* v2	2✔
55	density_weighted_velocity_3 = sqrt.(rho) .* v3	2✔
56
57	return compute_kinetic_energy_spectrum(density_weighted_velocity_1,	1✔
58	density_weighted_velocity_2,
59	density_weighted_velocity_3)
60	end
61
62	function interpolate_lgl_to_uniform_cartesian(u, mesh::TreeMesh{3},	1✔
63	equations::AbstractCompressibleEulerEquations,
64	solver::DGSEM, cache)
65	# Restrict to straightforward non AMR setups
66	leaf_cell_ids = leaf_cells(mesh.tree)	1✔
67	levels = mesh.tree.levels[leaf_cell_ids]	1✔
68	if !all(==(first(levels)), levels)	1✔
69	throw(ArgumentError("Non-uniform meshes are not supported yet"))	×
70	end
71
72	level = first(levels)	1✔
73	cells_per_dimension = 2^level	1✔
74
75	# Uses one uniform interpolation node per DGSEM solution node in each coordinate
76	# direction. A degree p element has p + 1 LGL nodes per direction, so this
77	# is the minimum Cartesian sampling matching the element-wise geometry
78	# but this can be increased to finer sampling sizes if needed
79	n_uniform_nodes = polydeg(solver) + 1	1✔
80	grid_points_per_dimension = n_uniform_nodes * cells_per_dimension	1✔
81
82	n_vars = nvariables(equations)	1✔
83	u_uniform = Array{eltype(u)}(undef, n_vars, grid_points_per_dimension,	1✔
84	grid_points_per_dimension,
85	grid_points_per_dimension)
86
87	# Interpolate from LGL nodes to cell-centered equidistant nodes in each element
88	dx_reference = 2 / n_uniform_nodes	1✔
89	nodes_out = collect(range(-1 + dx_reference / 2, 1 - dx_reference / 2,	1✔
90	length = n_uniform_nodes))
91	vandermonde = polynomial_interpolation_matrix(get_nodes(solver.basis), nodes_out)	1✔
92
93	# Element center coordinates determine where each interpolated block sits in the global grid
94	center = reshape(mesh.tree.center_level_0, :, 1)	1✔
95	normalized_coordinates = (mesh.tree.coordinates[:, leaf_cell_ids] .- center) ./	2✔
96	(mesh.tree.length_level_0 / 2)
97	dx_global = 2 / (n_uniform_nodes * cells_per_dimension)	1✔
98
99	for element in eachelement(solver, cache)	1✔
100	# Gather conservative nodal values on the reference LGL grid for the element
101	u_sample = get_node_vars(u, equations, solver, 1, 1, 1, element)	512✔
102	element_size = (n_vars, nnodes(solver), nnodes(solver),	512✔
103	nnodes(solver))
104	element_values = Array{eltype(u_sample)}(undef, element_size)	512✔
105	for k in eachnode(solver), j in eachnode(solver), i in eachnode(solver)	512✔
106	u_node = get_node_vars(u, equations, solver, i, j, k, element)	32,768✔
107	for variable in 1:n_vars	32,768✔
108	element_values[variable, i, j, k] = u_node[variable]	163,840✔
109	end	294,912✔
110	end	34,304✔
111	interpolated = multiply_dimensionwise(vandermonde, element_values)	512✔
112
113	# Each element is placed on the uniform grid assuming reference directions align with
114	# physical axes (ξ→x, η→y, ζ→z) and nodal values use the standard DGSEM tensor order along (ξ, η, ζ)
115	first_index = Vector{Int}(undef, 3)	512✔
116	for dim in 1:3	512✔
117	lower_left = normalized_coordinates[dim, element] -	1,536✔
118	(n_uniform_nodes - 1) / 2 * dx_global
119	first_index[dim] = round(Int,	1,536✔
120	(lower_left - (-1 + dx_global / 2)) / dx_global) +
121	1
122	end	2,560✔
123
124	# Writes the interpolated block onto the global grid for the larger output
125	# `r1`, `r2`, and `r3` are the global indices corresponding to `u_uniform` that this specific element's interpolated block fits within
126	# Essentially, `r1`, `r2`, and `r3` are the positions of the current element within `u_uniform`
127	# first_index[dim] refers to the first index of the nodes of this element within the array of global tensor nodes
128	# See the sketch in `spectral_analysis_2d.jl` for a sketch in 2d of how the local to global assembly works
129	r1 = first_index[1]:(first_index[1] + n_uniform_nodes - 1)	512✔
130	r2 = first_index[2]:(first_index[2] + n_uniform_nodes - 1)	512✔
131	r3 = first_index[3]:(first_index[3] + n_uniform_nodes - 1)	512✔
132	u_uniform[:, r1, r2, r3] .= interpolated[:, :, :, :]	1,024✔
133	end	1,023✔
134	return u_uniform	1✔
135	end
136
137	"""
138	compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{3}, equations,
139	dg::DGMultiSBP, cache)
140
141	Compute the energy spectrum for a 3D `DGMulti` finite-difference SBP solution whose
142	nodes already form a uniform Cartesian grid.
143	"""
144	function compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{3},	×
145	equations::AbstractCompressibleEulerEquations,
146	dg::DGMultiSBP, cache)
147	# Unpacks the primitive variables from the conservative state for FDSBP DGMulti solutions
148	u_values = parent(u)	×
149	n_points = length(u_values)	×
150	n = round(Int, n_points^(1 / 3))	×
151	q = cons2prim.(u_values, Ref(equations)) # q is the vector that contains the primitive variables for density and velocity converted from the conservative variables	×
152	rho = reshape(getindex.(q, 1), n, n, n)	×
153	density_weighted_velocity_1 = sqrt.(rho) .* reshape(getindex.(q, 2), n, n, n)	×
154	density_weighted_velocity_2 = sqrt.(rho) .* reshape(getindex.(q, 3), n, n, n)	×
155	density_weighted_velocity_3 = sqrt.(rho) .* reshape(getindex.(q, 4), n, n, n)	×
156
157	return compute_kinetic_energy_spectrum(density_weighted_velocity_1,	×
158	density_weighted_velocity_2,
159	density_weighted_velocity_3)
160	end
161	end # @muladd

trixi-framework / Trixi.jl / 30296365372

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