From 527f018b249d9b04f050ea738d1791f6c9d7aa49 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 8 Apr 2024 09:07:25 +0200 Subject: [PATCH 1/9] Use `eachindex` instead of `1:length` and compare against nothing with `===` --- .../src/files/non_periodic_boundaries.jl | 2 +- .../src/files/scalar_linear_advection_1d.jl | 4 ++-- .../elixir_advection_amr_solution_independent.jl | 2 +- .../elixir_advection_amr_solution_independent.jl | 2 +- .../elixir_advection_amr_solution_independent.jl | 2 +- ext/TrixiMakieExt.jl | 2 +- src/callbacks_step/amr.jl | 8 ++++---- .../euler_acoustics_coupling_dg2d.jl | 2 +- src/callbacks_step/time_series_dg_tree.jl | 16 ++++++++-------- .../time_series_dg_unstructured.jl | 12 ++++++------ src/meshes/abstract_tree.jl | 10 +++++----- src/meshes/parallel_tree_mesh.jl | 2 +- .../semidiscretization_coupled.jl | 4 ++-- .../semidiscretization_euler_acoustics.jl | 4 ++-- .../dgmulti/flux_differencing_gauss_sbp.jl | 4 ++-- src/solvers/dgsem_p4est/dg_parallel.jl | 6 +++--- src/solvers/dgsem_tree/dg_2d_parallel.jl | 2 +- src/solvers/dgsem_tree/dg_parallel.jl | 2 +- src/visualization/utilities.jl | 4 ++-- test/test_unit.jl | 2 +- utils/trixi2txt.jl | 8 ++++---- 21 files changed, 50 insertions(+), 50 deletions(-) diff --git a/docs/literate/src/files/non_periodic_boundaries.jl b/docs/literate/src/files/non_periodic_boundaries.jl index 3b238ad4533..8f0e320dfdc 100644 --- a/docs/literate/src/files/non_periodic_boundaries.jl +++ b/docs/literate/src/files/non_periodic_boundaries.jl @@ -99,7 +99,7 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), save_everystep=false, saveat=visnodes, callback=callbacks); using Plots -@gif for step in 1:length(sol.u) +@gif for step in eachindex(sol.u) plot(sol.u[step], semi, ylim=(-1.5, 1.5), legend=true, label="approximation", title="time t=$(round(sol.t[step], digits=5))") scatter!([0.0], [sum(boundary_condition(SVector(0.0), sol.t[step], equations))], label="boundary condition") end diff --git a/docs/literate/src/files/scalar_linear_advection_1d.jl b/docs/literate/src/files/scalar_linear_advection_1d.jl index 9b48f29d341..c7d55e26d2a 100644 --- a/docs/literate/src/files/scalar_linear_advection_1d.jl +++ b/docs/literate/src/files/scalar_linear_advection_1d.jl @@ -120,7 +120,7 @@ integral = sum(nodes.^3 .* weights) x = Matrix{Float64}(undef, length(nodes), n_elements) for element in 1:n_elements x_l = coordinates_min + (element - 1) * dx + dx/2 - for i in 1:length(nodes) + for i in eachindex(nodes) ξ = nodes[i] # nodes in [-1, 1] x[i, element] = x_l + dx/2 * ξ end @@ -417,7 +417,7 @@ dx = (coordinates_max - coordinates_min) / n_elements # length of one element x = Matrix{Float64}(undef, length(nodes), n_elements) for element in 1:n_elements x_l = -1 + (element - 1) * dx + dx/2 - for i in 1:length(nodes) # basis points in [-1, 1] + for i in eachindex(nodes) # basis points in [-1, 1] ξ = nodes[i] x[i, element] = x_l + dx/2 * ξ end diff --git a/examples/p4est_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/p4est_2d_dgsem/elixir_advection_amr_solution_independent.jl index 5a2537be4e6..a1ddc6a314f 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -33,7 +33,7 @@ function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, outer_distance = 1.85 #Iterate over all elements - for element in 1:length(alpha) + for element in eachindex(alpha) # Calculate periodic distance between cell and center. # This requires an uncurved mesh! coordinates = SVector(0.5 * (cache.elements.node_coordinates[1, 1, 1, element] + diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl index 0589e76a6a9..1ed08e1961b 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -32,7 +32,7 @@ function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, outer_distance = 1.85 # Iterate over all elements. - for element in 1:length(alpha) + for element in eachindex(alpha) # Calculate periodic distance between cell and center. # This requires an uncurved mesh! coordinates = SVector(0.5 * (cache.elements.node_coordinates[1, 1, 1, element] + diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl index 03a213689ec..c7412660b0c 100644 --- a/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -31,7 +31,7 @@ function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, outer_distance = 1.85 #Iterate over all elements - for element in 1:length(alpha) + for element in eachindex(alpha) #Calculate periodic distance between cell and center. cell_id = cache.elements.cell_ids[element] coordinates = mesh.tree.coordinates[1:2, cell_id] diff --git a/ext/TrixiMakieExt.jl b/ext/TrixiMakieExt.jl index 301a7656da9..1684ea2dc97 100644 --- a/ext/TrixiMakieExt.jl +++ b/ext/TrixiMakieExt.jl @@ -187,7 +187,7 @@ function iplot(pd::PlotData2DTriangulated; fig = Makie.Figure() # Set up options for the drop-down menu - menu_options = [zip(variable_names, 1:length(variable_names))...] + menu_options = [zip(variable_names, eachindex(variable_names))...] menu = Makie.Menu(fig, options = menu_options) # Initialize toggle switches for viewing the mesh diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 1ab65a3553e..45f03fba8fe 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -243,7 +243,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, @unpack to_refine, to_coarsen = amr_callback.amr_cache empty!(to_refine) empty!(to_coarsen) - for element in 1:length(lambda) + for element in eachindex(lambda) controller_value = lambda[element] if controller_value > 0 push!(to_refine, leaf_cell_ids[element]) @@ -307,7 +307,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, end # Extract only those parent cells for which all children should be coarsened - to_coarsen = collect(1:length(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] + to_coarsen = collect(eachindex(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] # Finally, coarsen mesh coarsened_original_cells = @trixi_timeit timer() "mesh" coarsen!(mesh.tree, @@ -395,7 +395,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, @unpack to_refine, to_coarsen = amr_callback.amr_cache empty!(to_refine) empty!(to_coarsen) - for element in 1:length(lambda) + for element in eachindex(lambda) controller_value = lambda[element] if controller_value > 0 push!(to_refine, leaf_cell_ids[element]) @@ -456,7 +456,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, end # Extract only those parent cells for which all children should be coarsened - to_coarsen = collect(1:length(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] + to_coarsen = collect(eachindex(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] # Finally, coarsen mesh coarsened_original_cells = @trixi_timeit timer() "mesh" coarsen!(mesh.tree, diff --git a/src/callbacks_step/euler_acoustics_coupling_dg2d.jl b/src/callbacks_step/euler_acoustics_coupling_dg2d.jl index 16fac4f2d8d..8a8bb893dcd 100644 --- a/src/callbacks_step/euler_acoustics_coupling_dg2d.jl +++ b/src/callbacks_step/euler_acoustics_coupling_dg2d.jl @@ -12,7 +12,7 @@ function calc_acoustic_sources!(acoustic_source_terms, u_euler, u_acoustics, dg::DGSEM, cache) acoustic_source_terms .= zero(eltype(acoustic_source_terms)) - @threaded for k in 1:length(coupled_element_ids) + @threaded for k in eachindex(coupled_element_ids) element = coupled_element_ids[k] for j in eachnode(dg), i in eachnode(dg) diff --git a/src/callbacks_step/time_series_dg_tree.jl b/src/callbacks_step/time_series_dg_tree.jl index 37d4e6ea705..0af1688a8ed 100644 --- a/src/callbacks_step/time_series_dg_tree.jl +++ b/src/callbacks_step/time_series_dg_tree.jl @@ -25,7 +25,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, mesh::TreeMesh, cell_id = cell_ids[element] # Iterate over coordinates - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Skip coordinates for which an element has already been found if element_ids[index] > 0 continue @@ -63,7 +63,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, wbary = barycentric_weights(nodes) - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Construct point x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh))) @@ -94,7 +94,7 @@ function record_state_at_points!(point_data, u, solution_variables, new_length = old_length + n_solution_variables # Loop over all points/elements that should be recorded - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Extract data array and element id data = point_data[index] element_id = element_ids[index] @@ -108,7 +108,7 @@ function record_state_at_points!(point_data, u, solution_variables, u_node = solution_variables(get_node_vars(u, equations, dg, i, element_id), equations) - for v in 1:length(u_node) + for v in eachindex(u_node) data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index]) end @@ -126,7 +126,7 @@ function record_state_at_points!(point_data, u, solution_variables, new_length = old_length + n_solution_variables # Loop over all points/elements that should be recorded - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Extract data array and element id data = point_data[index] element_id = element_ids[index] @@ -140,7 +140,7 @@ function record_state_at_points!(point_data, u, solution_variables, u_node = solution_variables(get_node_vars(u, equations, dg, i, j, element_id), equations) - for v in 1:length(u_node) + for v in eachindex(u_node) data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index] * interpolating_polynomials[j, 2, index]) @@ -159,7 +159,7 @@ function record_state_at_points!(point_data, u, solution_variables, new_length = old_length + n_solution_variables # Loop over all points/elements that should be recorded - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Extract data array and element id data = point_data[index] element_id = element_ids[index] @@ -173,7 +173,7 @@ function record_state_at_points!(point_data, u, solution_variables, u_node = solution_variables(get_node_vars(u, equations, dg, i, j, k, element_id), equations) - for v in 1:length(u_node) + for v in eachindex(u_node) data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index] * interpolating_polynomials[j, 2, index] diff --git a/src/callbacks_step/time_series_dg_unstructured.jl b/src/callbacks_step/time_series_dg_unstructured.jl index f6d1bb48f24..e9a18fcd3bb 100644 --- a/src/callbacks_step/time_series_dg_unstructured.jl +++ b/src/callbacks_step/time_series_dg_unstructured.jl @@ -31,7 +31,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, # Iterate over coordinates distances = zeros(eltype(mesh.corners), mesh.n_elements) indices = zeros(Int, mesh.n_elements, 2) - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Grab the current point for which the element needs found point = SVector(coordinates[1, index], coordinates[2, index]) @@ -77,7 +77,7 @@ function get_elements_by_coordinates!(element_ids, coordinates, # Loop through all the element candidates until we find a vector from the barycenter # to the surface that points in the same direction as the current `point` vector. # This then gives us the correct element. - for element in 1:length(candidates) + for element in eachindex(candidates) bary_center = SVector(bary_centers[1, candidates[element]], bary_centers[2, candidates[element]]) # Vector pointing from the barycenter toward the minimal `surface_point` @@ -118,7 +118,7 @@ function calc_minimum_surface_distance(point, node_coordinates, n = nnodes(dg) min_distance2 = Inf * ones(eltype(mesh.corners), length(mesh)) indices = zeros(Int, length(mesh), 2) - for k in 1:length(mesh) + for k in eachindex(mesh) # used to ensure that only boundary points are used on_surface = MVector(false, false) for j in 1:n @@ -153,7 +153,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, # Helper array for a straight-sided quadrilateral element corners = zeros(eltype(mesh.corners), 4, 2) - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Construct point x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh))) @@ -280,7 +280,7 @@ function record_state_at_points!(point_data, u, solution_variables, new_length = old_length + n_solution_variables # Loop over all points/elements that should be recorded - for index in 1:length(element_ids) + for index in eachindex(element_ids) # Extract data array and element id data = point_data[index] element_id = element_ids[index] @@ -294,7 +294,7 @@ function record_state_at_points!(point_data, u, solution_variables, u_node = solution_variables(get_node_vars(u, equations, dg, i, j, element_id), equations) - for v in 1:length(u_node) + for v in eachindex(u_node) data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index] * interpolating_polynomials[j, 2, index]) diff --git a/src/meshes/abstract_tree.jl b/src/meshes/abstract_tree.jl index 469189ff50c..7fe27b1966a 100644 --- a/src/meshes/abstract_tree.jl +++ b/src/meshes/abstract_tree.jl @@ -149,7 +149,7 @@ end function filter_leaf_cells(f, t::AbstractTree) filtered = Vector{Int}(undef, length(t)) count = 0 - for cell_id in 1:length(t) + for cell_id in eachindex(t) if is_leaf(t, cell_id) && f(cell_id) count += 1 filtered[count] = cell_id @@ -194,7 +194,7 @@ end # Store cell id in each cell to use for post-AMR analysis function reset_original_cell_ids!(t::AbstractTree) - t.original_cell_ids[1:length(t)] .= 1:length(t) + t.original_cell_ids[eachindex(t)] .= eachindex(t) end # Efficiently perform uniform refinement up to a given level (works only on mesh with a single cell) @@ -250,7 +250,7 @@ function init_neighbors!(t::AbstractTree, max_level = maximum_level(t)) # Initialize neighbors level by level for level in 1:max_level # Walk entire tree, starting from level 0 cell - for cell_id in 1:length(t) + for cell_id in eachindex(t) # Skip cells whose immediate children are already initialized *or* whose level is too high for this round if t.levels[cell_id] != level - 1 continue @@ -379,7 +379,7 @@ function refine!(t::AbstractTree, cell_ids, # Note: original_cell_ids contains the cell_id *before* refinement. At # refinement, the refined cell's original_cell_ids value has its sign flipped # to easily find it now. - refined_original_cells = @views(-t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0]) + refined_original_cells = @views(-t.original_cell_ids[eachindex(t)][t.original_cell_ids[eachindex(t)] .< 0]) # Check if count of refinement cells matches information in original_cell_ids @assert refinement_count==length(refined_original_cells) ("Mismatch in number of refined cells") @@ -605,7 +605,7 @@ function coarsen!(t::AbstractTree, cell_ids::AbstractArray{Int}) # Note: original_cell_ids contains the cell_id *before* coarsening. At # coarsening, the coarsened parent cell's original_cell_ids value has its sign flipped # to easily find it now. - @views coarsened_original_cells = (-t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0]) + @views coarsened_original_cells = (-t.original_cell_ids[eachindex(t)][t.original_cell_ids[eachindex(t)] .< 0]) # Check if count of coarsened cells matches information in original_cell_ids @assert n_coarsened==length(coarsened_original_cells) ("Mismatch in number of coarsened cells") diff --git a/src/meshes/parallel_tree_mesh.jl b/src/meshes/parallel_tree_mesh.jl index 050e419680c..da6ad665245 100644 --- a/src/meshes/parallel_tree_mesh.jl +++ b/src/meshes/parallel_tree_mesh.jl @@ -70,7 +70,7 @@ function partition!(mesh::ParallelTreeMesh; allow_coarsening = true) end end - @assert all(x -> x >= 0, mesh.tree.mpi_ranks[1:length(mesh.tree)]) + @assert all(x -> x >= 0, mesh.tree.mpi_ranks[eachindex(mesh.tree)]) @assert sum(mesh.n_cells_by_rank) == length(mesh.tree) return nothing diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index dc21dbe9a1e..9ed01644a52 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -32,14 +32,14 @@ function SemidiscretizationCoupled(semis...) # Number of coefficients for each semidiscretization n_coefficients = zeros(Int, length(semis)) - for i in 1:length(semis) + for i in eachindex(semis) _, equations, _, _ = mesh_equations_solver_cache(semis[i]) n_coefficients[i] = ndofs(semis[i]) * nvariables(equations) end # Compute range of coefficients associated with each semidiscretization and allocate coupled BCs u_indices = Vector{UnitRange{Int}}(undef, length(semis)) - for i in 1:length(semis) + for i in eachindex(semis) offset = sum(n_coefficients[1:(i - 1)]) + 1 u_indices[i] = range(offset, length = n_coefficients[i]) diff --git a/src/semidiscretization/semidiscretization_euler_acoustics.jl b/src/semidiscretization/semidiscretization_euler_acoustics.jl index 173523ff892..286315fb960 100644 --- a/src/semidiscretization/semidiscretization_euler_acoustics.jl +++ b/src/semidiscretization/semidiscretization_euler_acoustics.jl @@ -103,7 +103,7 @@ function precompute_weights(source_region, weights, coupled_element_ids, equatio (nnodes(dg), nnodes(dg), length(coupled_element_ids))) - @threaded for k in 1:length(coupled_element_ids) + @threaded for k in eachindex(coupled_element_ids) element = coupled_element_ids[k] for j in eachnode(dg), i in eachnode(dg) x = get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, @@ -197,7 +197,7 @@ function add_acoustic_source_terms!(du_acoustics, acoustic_source_terms, source_ coupled_element_ids, mesh::TreeMesh{2}, equations, dg::DGSEM, cache) - @threaded for k in 1:length(coupled_element_ids) + @threaded for k in eachindex(coupled_element_ids) element = coupled_element_ids[k] for j in eachnode(dg), i in eachnode(dg) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index 9059caf87f6..cab02f9e565 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -63,7 +63,7 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, # Permutation of indices in a tensor product form num_faces = StartUpDG.num_faces(rd.element_type) - indices = reshape(1:length(rd.rf), nnodes_1d, num_faces) + indices = reshape(eachindex(rd.rf), nnodes_1d, num_faces) face_indices_tensor_product = zeros(Int, 2, nnodes_1d, ndims(rd.element_type)) for i in 1:nnodes_1d # loop over nodes in one face face_indices_tensor_product[:, i, 1] .= indices[i, 1:2] @@ -93,7 +93,7 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, # Permutation of indices in a tensor product form num_faces = StartUpDG.num_faces(rd.element_type) - indices = reshape(1:length(rd.rf), nnodes_1d, nnodes_1d, num_faces) + indices = reshape(eachindex(rd.rf), nnodes_1d, nnodes_1d, num_faces) face_indices_tensor_product = zeros(Int, 2, nnodes_1d, nnodes_1d, ndims(rd.element_type)) for j in 1:nnodes_1d, i in 1:nnodes_1d # loop over nodes in one face diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index eaa6ab5cee2..43649a5dc6e 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -49,7 +49,7 @@ function start_mpi_send!(mpi_cache::P4estMPICache, mesh, equations, dg, cache) data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1) n_small_elements = 2^(ndims(mesh) - 1) - for d in 1:length(mpi_cache.mpi_neighbor_ranks) + for d in eachindex(mpi_cache.mpi_neighbor_ranks) send_buffer = mpi_cache.mpi_send_buffers[d] for (index, interface) in enumerate(mpi_cache.mpi_neighbor_interfaces[d]) @@ -505,7 +505,7 @@ function exchange_normal_directions!(mpi_mortars, mpi_cache, # Create buffers and requests send_buffers = Vector{Vector{RealT}}(undef, length(mpi_neighbor_mortars)) recv_buffers = Vector{Vector{RealT}}(undef, length(mpi_neighbor_mortars)) - for index in 1:length(mpi_neighbor_mortars) + for index in eachindex(mpi_neighbor_mortars) send_buffers[index] = Vector{RealT}(undef, length(mpi_neighbor_mortars[index]) * n_small_elements * data_size) @@ -519,7 +519,7 @@ function exchange_normal_directions!(mpi_mortars, mpi_cache, recv_requests = Vector{MPI.Request}(undef, length(mpi_neighbor_mortars)) # Fill send buffers - for d in 1:length(mpi_neighbor_ranks) + for d in eachindex(mpi_neighbor_ranks) send_buffer = send_buffers[d] for (index, mortar) in enumerate(mpi_neighbor_mortars[d]) diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index 157d462aa2f..746a9802d66 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -60,7 +60,7 @@ end function start_mpi_send!(mpi_cache::MPICache, mesh, equations, dg, cache) data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1) - for d in 1:length(mpi_cache.mpi_neighbor_ranks) + for d in eachindex(mpi_cache.mpi_neighbor_ranks) send_buffer = mpi_cache.mpi_send_buffers[d] for (index, interface) in enumerate(mpi_cache.mpi_neighbor_interfaces[d]) diff --git a/src/solvers/dgsem_tree/dg_parallel.jl b/src/solvers/dgsem_tree/dg_parallel.jl index c614fe0d0e6..f1e78b2b64a 100644 --- a/src/solvers/dgsem_tree/dg_parallel.jl +++ b/src/solvers/dgsem_tree/dg_parallel.jl @@ -13,7 +13,7 @@ function init_mpi_data_structures(mpi_neighbor_interfaces, mpi_neighbor_mortars, n_small_elements = 2^(n_dims - 1) mpi_send_buffers = Vector{Vector{uEltype}}(undef, length(mpi_neighbor_interfaces)) mpi_recv_buffers = Vector{Vector{uEltype}}(undef, length(mpi_neighbor_interfaces)) - for index in 1:length(mpi_neighbor_interfaces) + for index in eachindex(mpi_neighbor_interfaces) mpi_send_buffers[index] = Vector{uEltype}(undef, length(mpi_neighbor_interfaces[index]) * data_size + diff --git a/src/visualization/utilities.jl b/src/visualization/utilities.jl index 05457395ac0..c1108128c92 100644 --- a/src/visualization/utilities.jl +++ b/src/visualization/utilities.jl @@ -495,7 +495,7 @@ function cell2node(cell_centered_data) resolution_in, _ = size(first(cell_centered_data)) resolution_out = resolution_in + 1 node_centered_data = [Matrix{Float64}(undef, resolution_out, resolution_out) - for _ in 1:length(cell_centered_data)] + for _ in eachindex(cell_centered_data)] for (cell_data, node_data) in zip(cell_centered_data, node_centered_data) # Fill center with original data @@ -1545,7 +1545,7 @@ end # Create an axis. function axis_curve(nodes_x, nodes_y, nodes_z, slice, point, n_points) - if n_points == nothing + if n_points === nothing n_points = 64 end dimensions = length(point) diff --git a/test/test_unit.jl b/test/test_unit.jl index 79950f58d59..9a4ae2ede13 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -284,7 +284,7 @@ end function MyContainer(data, capacity) c = MyContainer(Vector{Int}(undef, capacity + 1), capacity, length(data), capacity + 1) - c.data[1:length(data)] .= data + c.data[eachindex(data)] .= data return c end MyContainer(data::AbstractArray) = MyContainer(data, length(data)) diff --git a/utils/trixi2txt.jl b/utils/trixi2txt.jl index 12a3d46760e..52ee904d2f6 100644 --- a/utils/trixi2txt.jl +++ b/utils/trixi2txt.jl @@ -86,7 +86,7 @@ function trixi2txt(filename::AbstractString...; "maximum supported level $max_supported_level") end max_available_nodes_per_finest_element = 2^(max_supported_level - max_level) - if nvisnodes == nothing + if nvisnodes === nothing max_nvisnodes = 2 * n_nodes elseif nvisnodes == 0 max_nvisnodes = n_nodes @@ -137,9 +137,9 @@ function trixi2txt(filename::AbstractString...; println(io) # Data - for idx in 1:length(xs) + for idx in eachindex(xs) @printf(io, "%+10.8e", xs[idx]) - for variable_id in 1:length(variables) + for variable_id in eachindex(variables) @printf(io, " %+10.8e ", node_centered_data[idx, variable_id]) end println(io) @@ -199,7 +199,7 @@ function read_meshfile(filename::String) # Extract leaf cells (= cells to be plotted) and contract all other arrays accordingly leaf_cells = similar(levels) n_cells = 0 - for cell_id in 1:length(levels) + for cell_id in eachindex(levels) if sum(child_ids[:, cell_id]) > 0 continue end From 2024dae1f17073941e52c787c31337619deef256 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 8 Apr 2024 09:28:24 +0200 Subject: [PATCH 2/9] revert some places --- src/meshes/abstract_tree.jl | 10 +++++----- src/meshes/parallel_tree_mesh.jl | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/meshes/abstract_tree.jl b/src/meshes/abstract_tree.jl index 7fe27b1966a..469189ff50c 100644 --- a/src/meshes/abstract_tree.jl +++ b/src/meshes/abstract_tree.jl @@ -149,7 +149,7 @@ end function filter_leaf_cells(f, t::AbstractTree) filtered = Vector{Int}(undef, length(t)) count = 0 - for cell_id in eachindex(t) + for cell_id in 1:length(t) if is_leaf(t, cell_id) && f(cell_id) count += 1 filtered[count] = cell_id @@ -194,7 +194,7 @@ end # Store cell id in each cell to use for post-AMR analysis function reset_original_cell_ids!(t::AbstractTree) - t.original_cell_ids[eachindex(t)] .= eachindex(t) + t.original_cell_ids[1:length(t)] .= 1:length(t) end # Efficiently perform uniform refinement up to a given level (works only on mesh with a single cell) @@ -250,7 +250,7 @@ function init_neighbors!(t::AbstractTree, max_level = maximum_level(t)) # Initialize neighbors level by level for level in 1:max_level # Walk entire tree, starting from level 0 cell - for cell_id in eachindex(t) + for cell_id in 1:length(t) # Skip cells whose immediate children are already initialized *or* whose level is too high for this round if t.levels[cell_id] != level - 1 continue @@ -379,7 +379,7 @@ function refine!(t::AbstractTree, cell_ids, # Note: original_cell_ids contains the cell_id *before* refinement. At # refinement, the refined cell's original_cell_ids value has its sign flipped # to easily find it now. - refined_original_cells = @views(-t.original_cell_ids[eachindex(t)][t.original_cell_ids[eachindex(t)] .< 0]) + refined_original_cells = @views(-t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0]) # Check if count of refinement cells matches information in original_cell_ids @assert refinement_count==length(refined_original_cells) ("Mismatch in number of refined cells") @@ -605,7 +605,7 @@ function coarsen!(t::AbstractTree, cell_ids::AbstractArray{Int}) # Note: original_cell_ids contains the cell_id *before* coarsening. At # coarsening, the coarsened parent cell's original_cell_ids value has its sign flipped # to easily find it now. - @views coarsened_original_cells = (-t.original_cell_ids[eachindex(t)][t.original_cell_ids[eachindex(t)] .< 0]) + @views coarsened_original_cells = (-t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0]) # Check if count of coarsened cells matches information in original_cell_ids @assert n_coarsened==length(coarsened_original_cells) ("Mismatch in number of coarsened cells") diff --git a/src/meshes/parallel_tree_mesh.jl b/src/meshes/parallel_tree_mesh.jl index da6ad665245..050e419680c 100644 --- a/src/meshes/parallel_tree_mesh.jl +++ b/src/meshes/parallel_tree_mesh.jl @@ -70,7 +70,7 @@ function partition!(mesh::ParallelTreeMesh; allow_coarsening = true) end end - @assert all(x -> x >= 0, mesh.tree.mpi_ranks[eachindex(mesh.tree)]) + @assert all(x -> x >= 0, mesh.tree.mpi_ranks[1:length(mesh.tree)]) @assert sum(mesh.n_cells_by_rank) == length(mesh.tree) return nothing From 880a642c8c85a5ff3afeaf3451916a11525d2421 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 8 Apr 2024 09:30:24 +0200 Subject: [PATCH 3/9] revert for solvers --- .../dgmulti/flux_differencing_gauss_sbp.jl | 1162 ++++++++++------- src/solvers/dgsem_p4est/dg_parallel.jl | 6 +- src/solvers/dgsem_tree/dg_parallel.jl | 2 +- 3 files changed, 678 insertions(+), 492 deletions(-) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index cab02f9e565..469189ff50c 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -1,593 +1,779 @@ +# 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 -# ========= GaussSBP approximation types ============ -# Note: we define type aliases outside of the @muladd block to avoid Revise breaking when code -# inside the @muladd block is edited. See https://github.com/trixi-framework/Trixi.jl/issues/801 -# for more details. +abstract type AbstractTree{NDIMS} <: AbstractContainer end -# `GaussSBP` is a type alias for a StartUpDG type (e.g., Gauss nodes on quads/hexes) -const GaussSBP = Polynomial{Gauss} +# Type traits to obtain dimension +@inline Base.ndims(::AbstractTree{NDIMS}) where {NDIMS} = NDIMS -function tensor_product_quadrature(element_type::Line, r1D, w1D) - return r1D, w1D +# Auxiliary methods to allow semantic queries on the tree +# Check whether cell has parent cell +has_parent(t::AbstractTree, cell_id::Int) = t.parent_ids[cell_id] > 0 + +# Count number of children for a given cell +function n_children(t::AbstractTree, cell_id::Int) + count(x -> (x > 0), @view t.child_ids[:, cell_id]) end -function tensor_product_quadrature(element_type::Quad, r1D, w1D) - sq, rq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D)) - ws, wr = vec.(StartUpDG.NodesAndModes.meshgrid(w1D)) - wq = wr .* ws - return rq, sq, wq +# Check whether cell has any child cell +has_children(t::AbstractTree, cell_id::Int) = n_children(t, cell_id) > 0 + +# Check whether cell is leaf cell +is_leaf(t::AbstractTree, cell_id::Int) = !has_children(t, cell_id) + +# Check whether cell has specific child cell +has_child(t::AbstractTree, cell_id::Int, child::Int) = t.child_ids[child, cell_id] > 0 + +# Check if cell has a neighbor at the same refinement level in the given direction +function has_neighbor(t::AbstractTree, cell_id::Int, direction::Int) + t.neighbor_ids[direction, cell_id] > 0 end -function tensor_product_quadrature(element_type::Hex, r1D, w1D) - rq, sq, tq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D, r1D, r1D)) - wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(w1D, w1D, w1D)) - wq = wr .* ws .* wt - return rq, sq, tq, wq +# Check if cell has a coarse neighbor, i.e., with one refinement level lower +function has_coarse_neighbor(t::AbstractTree, cell_id::Int, direction::Int) + return has_parent(t, cell_id) && has_neighbor(t, t.parent_ids[cell_id], direction) end -# type parameters for `TensorProductFaceOperator`. -abstract type AbstractGaussOperator end -struct Interpolation <: AbstractGaussOperator end -# - `Projection{ScaleByFaceWeights=Static.False()}` corresponds to the operator `projection_matrix_gauss_to_face = M \ Vf'`, -# which is used in `VolumeIntegralFluxDifferencing`. -# - `Projection{ScaleByFaceWeights=Static.True()}` corresponds to the quadrature-based lifting -# operator `LIFT = M \ (Vf' * diagm(rd.wf))`, which is used in `SurfaceIntegralWeakForm` -struct Projection{ScaleByFaceWeights} <: AbstractGaussOperator end +# Check if cell has any neighbor (same-level or lower-level) +function has_any_neighbor(t::AbstractTree, cell_id::Int, direction::Int) + return has_neighbor(t, cell_id, direction) || + has_coarse_neighbor(t, cell_id, direction) +end -# used to dispatch for different Gauss interpolation operators -abstract type AbstractTensorProductGaussOperator end +# Check if cell is own cell, i.e., belongs to this MPI rank +is_own_cell(t::AbstractTree, cell_id) = true -# TensorProductGaussFaceOperator{Tmat, Ti} -# -# Data for performing tensor product interpolation from volume nodes to face nodes. -struct TensorProductGaussFaceOperator{NDIMS, OperatorType <: AbstractGaussOperator, - Tmat, Tweights, Tfweights, Tindices} <: - AbstractTensorProductGaussOperator - interp_matrix_gauss_to_face_1d::Tmat - inv_volume_weights_1d::Tweights - face_weights::Tfweights - face_indices_tensor_product::Tindices - nnodes_1d::Int - nfaces::Int -end +# Return cell length for a given level +length_at_level(t::AbstractTree, level::Int) = t.length_level_0 / 2^level -# constructor for a 2D operator -function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, - dg::DGMulti{2, Quad, GaussSBP}) - rd = dg.basis +# Return cell length for a given cell +length_at_cell(t::AbstractTree, cell_id::Int) = length_at_level(t, t.levels[cell_id]) - rq1D, wq1D = StartUpDG.gauss_quad(0, 0, polydeg(dg)) - interp_matrix_gauss_to_face_1d = polynomial_interpolation_matrix(rq1D, [-1; 1]) +# Return minimum level of any leaf cell +minimum_level(t::AbstractTree) = minimum(t.levels[leaf_cells(t)]) - nnodes_1d = length(rq1D) +# Return maximum level of any leaf cell +maximum_level(t::AbstractTree) = maximum(t.levels[leaf_cells(t)]) - # Permutation of indices in a tensor product form - num_faces = StartUpDG.num_faces(rd.element_type) - indices = reshape(eachindex(rd.rf), nnodes_1d, num_faces) - face_indices_tensor_product = zeros(Int, 2, nnodes_1d, ndims(rd.element_type)) - for i in 1:nnodes_1d # loop over nodes in one face - face_indices_tensor_product[:, i, 1] .= indices[i, 1:2] - face_indices_tensor_product[:, i, 2] .= indices[i, 3:4] - end +# Check if tree is periodic +isperiodic(t::AbstractTree) = all(t.periodicity) +isperiodic(t::AbstractTree, dimension) = t.periodicity[dimension] + +# Auxiliary methods for often-required calculations +# Number of potential child cells +n_children_per_cell(::AbstractTree{NDIMS}) where {NDIMS} = 2^NDIMS - T_op = typeof(operator) - Tm = typeof(interp_matrix_gauss_to_face_1d) - Tw = typeof(inv.(wq1D)) - Tf = typeof(rd.wf) - Ti = typeof(face_indices_tensor_product) - return TensorProductGaussFaceOperator{2, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d, - inv.(wq1D), rd.wf, - face_indices_tensor_product, - nnodes_1d, num_faces) +# Number of directions +# +# Directions are indicated by numbers from 1 to 2*ndims: +# 1 -> -x +# 2 -> +x +# 3 -> -y +# 4 -> +y +# 5 -> -z +# 6 -> +z +@inline n_directions(::AbstractTree{NDIMS}) where {NDIMS} = 2 * NDIMS +# TODO: Taal performance, 1:n_directions(tree) vs. Base.OneTo(n_directions(tree)) vs. SOneTo(n_directions(tree)) +""" + eachdirection(tree::AbstractTree) + +Return an iterator over the indices that specify the location in relevant data structures +for the directions in `AbstractTree`. +In particular, not the directions themselves are returned. +""" +@inline eachdirection(tree::AbstractTree) = Base.OneTo(n_directions(tree)) + +# For a given direction, return its opposite direction +# +# dir -> opp +# 1 -> 2 +# 2 -> 1 +# 3 -> 4 +# 4 -> 3 +# 5 -> 6 +# 6 -> 5 +opposite_direction(direction::Int) = direction + 1 - 2 * ((direction + 1) % 2) + +# For a given child position (from 1 to 8) and dimension (from 1 to 3), +# calculate a child cell's position relative to its parent cell. +# +# Essentially calculates the following +# dim=1 dim=2 dim=3 +# child x y z +# 1 - - - +# 2 + - - +# 3 - + - +# 4 + + - +# 5 - - + +# 6 + - + +# 7 - + + +# 8 + + + +# child_sign(child::Int, dim::Int) = 1 - 2 * (div(child + 2^(dim - 1) - 1, 2^(dim-1)) % 2) +# Since we use only a fixed number of dimensions, we use a lookup table for improved performance. +const _child_signs = [-1 -1 -1; + +1 -1 -1; + -1 +1 -1; + +1 +1 -1; + -1 -1 +1; + +1 -1 +1; + -1 +1 +1; + +1 +1 +1] +child_sign(child::Int, dim::Int) = _child_signs[child, dim] + +# For each child position (1 to 8) and a given direction (from 1 to 6), return +# neighboring child position. +const _adjacent_child_ids = [2 2 3 3 5 5; + 1 1 4 4 6 6; + 4 4 1 1 7 7; + 3 3 2 2 8 8; + 6 6 7 7 1 1; + 5 5 8 8 2 2; + 8 8 5 5 3 3; + 7 7 6 6 4 4] +adjacent_child(child::Int, direction::Int) = _adjacent_child_ids[child, direction] + +# For each child position (1 to 8) and a given direction (from 1 to 6), return +# if neighbor is a sibling +function has_sibling(child::Int, direction::Int) + return (child_sign(child, div(direction + 1, 2)) * (-1)^(direction - 1)) > 0 end -# constructor for a 3D operator -function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, - dg::DGMulti{3, Hex, GaussSBP}) - rd = dg.basis - - rq1D, wq1D = StartUpDG.gauss_quad(0, 0, polydeg(dg)) - interp_matrix_gauss_to_face_1d = polynomial_interpolation_matrix(rq1D, [-1; 1]) - - nnodes_1d = length(rq1D) - - # Permutation of indices in a tensor product form - num_faces = StartUpDG.num_faces(rd.element_type) - indices = reshape(eachindex(rd.rf), nnodes_1d, nnodes_1d, num_faces) - face_indices_tensor_product = zeros(Int, 2, nnodes_1d, nnodes_1d, - ndims(rd.element_type)) - for j in 1:nnodes_1d, i in 1:nnodes_1d # loop over nodes in one face - face_indices_tensor_product[:, i, j, 1] .= indices[i, j, 1:2] - face_indices_tensor_product[:, i, j, 2] .= indices[i, j, 3:4] - face_indices_tensor_product[:, i, j, 3] .= indices[i, j, 5:6] +# Obtain leaf cells that fulfill a given criterion. +# +# The function `f` is passed the cell id of each leaf cell +# as an argument. +function filter_leaf_cells(f, t::AbstractTree) + filtered = Vector{Int}(undef, length(t)) + count = 0 + for cell_id in 1:length(t) + if is_leaf(t, cell_id) && f(cell_id) + count += 1 + filtered[count] = cell_id + end end - T_op = typeof(operator) - Tm = typeof(interp_matrix_gauss_to_face_1d) - Tw = typeof(inv.(wq1D)) - Tf = typeof(rd.wf) - Ti = typeof(face_indices_tensor_product) - return TensorProductGaussFaceOperator{3, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d, - inv.(wq1D), rd.wf, - face_indices_tensor_product, - nnodes_1d, num_faces) + return filtered[1:count] end -# specialize behavior of `mul_by!(A)` where `A isa TensorProductGaussFaceOperator)` -@inline function mul_by!(A::AbstractTensorProductGaussOperator) - return (out, x) -> tensor_product_gauss_face_operator!(out, A, x) +# Return an array with the ids of all leaf cells +leaf_cells(t::AbstractTree) = filter_leaf_cells((cell_id) -> true, t) + +# Return an array with the ids of all leaf cells for a given rank +leaf_cells_by_rank(t::AbstractTree, rank) = leaf_cells(t) + +# Return an array with the ids of all local leaf cells +local_leaf_cells(t::AbstractTree) = leaf_cells(t) + +# Count the number of leaf cells. +count_leaf_cells(t::AbstractTree) = length(leaf_cells(t)) + +@inline function cell_coordinates(t::AbstractTree{NDIMS}, cell) where {NDIMS} + SVector(ntuple(d -> t.coordinates[d, cell], Val(NDIMS))) end -@inline function tensor_product_gauss_face_operator!(out::AbstractMatrix, - A::AbstractTensorProductGaussOperator, - x::AbstractMatrix) - @threaded for col in Base.OneTo(size(out, 2)) - tensor_product_gauss_face_operator!(view(out, :, col), A, view(x, :, col)) +@inline function set_cell_coordinates!(t::AbstractTree{NDIMS}, coords, + cell) where {NDIMS} + for d in 1:NDIMS + t.coordinates[d, cell] = coords[d] end end -# 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 +# Determine if point is located inside cell +function is_point_in_cell(t::AbstractTree, point_coordinates, cell_id) + cell_length = length_at_cell(t, cell_id) + cell_coordinates_ = cell_coordinates(t, cell_id) + min_coordinates = cell_coordinates_ .- cell_length / 2 + max_coordinates = cell_coordinates_ .+ cell_length / 2 -#! format: off -# Interpolates values from volume Gauss nodes to face nodes on one element. -@inline function tensor_product_gauss_face_operator!(out::AbstractVector, - A::TensorProductGaussFaceOperator{2, Interpolation}, - x_in::AbstractVector) -#! format: on - (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A - (; nnodes_1d) = A - - fill!(out, zero(eltype(out))) - - # for 2D GaussSBP nodes, the indexing is first in x, then in y - x = reshape(x_in, nnodes_1d, nnodes_1d) - - # interpolation in the x-direction - @turbo for i in Base.OneTo(nnodes_1d) # loop over nodes in a face - index_left = face_indices_tensor_product[1, i, 1] - index_right = face_indices_tensor_product[2, i, 1] - for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes - out[index_left] = out[index_left] + - interp_matrix_gauss_to_face_1d[1, jj] * x[jj, i] - out[index_right] = out[index_right] + - interp_matrix_gauss_to_face_1d[2, jj] * x[jj, i] - end - end + return all(min_coordinates .<= point_coordinates .<= max_coordinates) +end - # interpolation in the y-direction - @turbo for i in Base.OneTo(nnodes_1d) # loop over nodes in a face - index_left = face_indices_tensor_product[1, i, 2] - index_right = face_indices_tensor_product[2, i, 2] - for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes - out[index_left] = out[index_left] + - interp_matrix_gauss_to_face_1d[1, jj] * x[i, jj] - out[index_right] = out[index_right] + - interp_matrix_gauss_to_face_1d[2, jj] * x[i, jj] - end - end +# Store cell id in each cell to use for post-AMR analysis +function reset_original_cell_ids!(t::AbstractTree) + t.original_cell_ids[1:length(t)] .= 1:length(t) end -# Interpolates values from volume Gauss nodes to face nodes on one element. -#! format: off -@inline function tensor_product_gauss_face_operator!(out::AbstractVector, - A::TensorProductGaussFaceOperator{3, Interpolation}, - x::AbstractVector) -#! format: on - (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A - (; nnodes_1d) = A - - fill!(out, zero(eltype(out))) - - # for 3D GaussSBP nodes, the indexing is first in y, then x, then z. - x = reshape(x, nnodes_1d, nnodes_1d, nnodes_1d) - - # interpolation in the y-direction - @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face - index_left = face_indices_tensor_product[1, i, j, 2] - index_right = face_indices_tensor_product[2, i, j, 2] - for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes - out[index_left] = out[index_left] + - interp_matrix_gauss_to_face_1d[1, jj] * x[jj, i, j] - out[index_right] = out[index_right] + - interp_matrix_gauss_to_face_1d[2, jj] * x[jj, i, j] - end - end +# Efficiently perform uniform refinement up to a given level (works only on mesh with a single cell) +function refine_uniformly!(t::AbstractTree, max_level) + @assert length(t)==1 "efficient uniform refinement only works for a newly created tree" + @assert max_level>=0 "the uniform refinement level must be non-zero" - # interpolation in the x-direction - @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face - index_left = face_indices_tensor_product[1, i, j, 1] - index_right = face_indices_tensor_product[2, i, j, 1] - for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes - out[index_left] = out[index_left] + - interp_matrix_gauss_to_face_1d[1, jj] * x[i, jj, j] - out[index_right] = out[index_right] + - interp_matrix_gauss_to_face_1d[2, jj] * x[i, jj, j] - end + # Calculate size of final tree and resize tree + total_length = 1 + for level in 1:max_level + total_length += n_children_per_cell(t)^level end + resize!(t, total_length) - # interpolation in the z-direction - @turbo for i in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d) # loop over nodes in a face - index_left = face_indices_tensor_product[1, i, j, 3] - index_right = face_indices_tensor_product[2, i, j, 3] - for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes - # The ordering (i,j) -> (j,i) needs to be reversed for this last face. - # This is due to way we define face nodes for Hex() types in StartUpDG.jl. - out[index_left] = out[index_left] + - interp_matrix_gauss_to_face_1d[1, jj] * x[j, i, jj] - out[index_right] = out[index_right] + - interp_matrix_gauss_to_face_1d[2, jj] * x[j, i, jj] - end - end + # Traverse tree to set parent-child relationships + init_children!(t, 1, max_level) + + # Set all neighbor relationships + init_neighbors!(t, max_level) end -# Projects face node values to volume Gauss nodes on one element. -#! format: off -@inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector, - A::TensorProductGaussFaceOperator{2, Projection{ApplyFaceWeights}}, - x::AbstractVector) where {ApplyFaceWeights} -#! format: on - (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A - (; inv_volume_weights_1d, nnodes_1d) = A - - fill!(out_vec, zero(eltype(out_vec))) - - # As of Julia 1.9, Base.ReshapedArray does not produce allocations when setting values. - # Thus, Base.ReshapedArray should be used if you are setting values in the array. - # `reshape` is fine if you are only accessing values. - # Note that, for 2D GaussSBP nodes, the indexing is first in x, then y - out = Base.ReshapedArray(out_vec, (nnodes_1d, nnodes_1d), ()) - - if ApplyFaceWeights == true - @turbo for i in eachindex(x) - x[i] = x[i] * A.face_weights[i] +# Recursively initialize children up to level `max_level` in depth-first ordering, starting with +# cell `cell_id` and set all information except neighbor relations (see `init_neighbors!`). +# +# Return the number of offspring of the initialized cell plus one +function init_children!(t::AbstractTree, cell_id, max_level) + # Stop recursion if max_level has been reached + if t.levels[cell_id] >= max_level + return 1 + else + # Initialize each child cell, counting the total number of offspring + n_offspring = 0 + for child in 1:n_children_per_cell(t) + # Get cell id of child + child_id = cell_id + 1 + n_offspring + + # Initialize child cell (except neighbors) + init_child!(t, cell_id, child, child_id) + + # Recursively initialize child cell + n_offspring += init_children!(t, child_id, max_level) end - end - # interpolation in the x-direction - @turbo for i in Base.OneTo(nnodes_1d) # loop over face nodes - index_left = face_indices_tensor_product[1, i, 1] - index_right = face_indices_tensor_product[2, i, 1] - for jj in Base.OneTo(nnodes_1d) # loop over a line of volume nodes - out[jj, i] = out[jj, i] + - interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] - out[jj, i] = out[jj, i] + - interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] - end + return n_offspring + 1 end +end - # interpolation in the y-direction - @turbo for i in Base.OneTo(nnodes_1d) - index_left = face_indices_tensor_product[1, i, 2] - index_right = face_indices_tensor_product[2, i, 2] - # loop over a line of volume nodes - for jj in Base.OneTo(nnodes_1d) - out[i, jj] = out[i, jj] + - interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] - out[i, jj] = out[i, jj] + - interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] +# Iteratively set all neighbor relations, starting at an initialized level 0 cell. Assume that +# parent-child relations have already been initialized (see `init_children!`). +function init_neighbors!(t::AbstractTree, max_level = maximum_level(t)) + @assert all(n >= 0 for n in t.neighbor_ids[:, 1]) "level 0 cell neighbors must be initialized" + + # Initialize neighbors level by level + for level in 1:max_level + # Walk entire tree, starting from level 0 cell + for cell_id in 1:length(t) + # Skip cells whose immediate children are already initialized *or* whose level is too high for this round + if t.levels[cell_id] != level - 1 + continue + end + + # Iterate over children and set neighbor information + for child in 1:n_children_per_cell(t) + child_id = t.child_ids[child, cell_id] + init_child_neighbors!(t, cell_id, child, child_id) + end end end - # apply inv(M) - @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) - out[i, j] = out[i, j] * inv_volume_weights_1d[i] * inv_volume_weights_1d[j] - end + return nothing end -# Interpolates values from volume Gauss nodes to face nodes on one element. -#! format: off -@inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector, - A::TensorProductGaussFaceOperator{3, Projection{ApplyFaceWeights}}, - x::AbstractVector) where {ApplyFaceWeights} -#! format: on - @unpack interp_matrix_gauss_to_face_1d, face_indices_tensor_product = A - @unpack inv_volume_weights_1d, nnodes_1d, nfaces = A - - fill!(out_vec, zero(eltype(out_vec))) - - # As of Julia 1.9, Base.ReshapedArray does not produce allocations when setting values. - # Thus, Base.ReshapedArray should be used if you are setting values in the array. - # `reshape` is fine if you are only accessing values. - # Note that, for 3D GaussSBP nodes, the indexing is first in y, then x, then z. - out = Base.ReshapedArray(out_vec, (nnodes_1d, nnodes_1d, nnodes_1d), ()) - - if ApplyFaceWeights == true - @turbo for i in eachindex(x) - x[i] = x[i] * A.face_weights[i] +# Initialize the neighbors of child cell `child_id` based on parent cell `cell_id` +function init_child_neighbors!(t::AbstractTree, cell_id, child, child_id) + t.neighbor_ids[:, child_id] .= zero(eltype(t.neighbor_ids)) + for direction in eachdirection(t) + # If neighbor is a sibling, establish one-sided connectivity + # Note: two-sided is not necessary, as each sibling will do this + if has_sibling(child, direction) + adjacent = adjacent_child(child, direction) + neighbor_id = t.child_ids[adjacent, cell_id] + + t.neighbor_ids[direction, child_id] = neighbor_id + continue end - end - # interpolation in the y-direction - @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face - index_left = face_indices_tensor_product[1, i, j, 2] - index_right = face_indices_tensor_product[2, i, j, 2] - for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes - out[jj, i, j] = out[jj, i, j] + - interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] - out[jj, i, j] = out[jj, i, j] + - interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] + # Skip if original cell does have no neighbor in direction + if !has_neighbor(t, cell_id, direction) + continue end - end - # interpolation in the x-direction - @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face - index_left = face_indices_tensor_product[1, i, j, 1] - index_right = face_indices_tensor_product[2, i, j, 1] - for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes - out[i, jj, j] = out[i, jj, j] + - interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] - out[i, jj, j] = out[i, jj, j] + - interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] + # Otherwise, check if neighbor has children - if not, skip again + neighbor_id = t.neighbor_ids[direction, cell_id] + if !has_children(t, neighbor_id) + continue end - end - # interpolation in the z-direction - @turbo for i in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d) # loop over nodes in a face - index_left = face_indices_tensor_product[1, i, j, 3] - index_right = face_indices_tensor_product[2, i, j, 3] - for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes - # The ordering (i,j) -> (j,i) needs to be reversed for this last face. - # This is due to way we define face nodes for Hex() types in StartUpDG.jl. - out[j, i, jj] = out[j, i, jj] + - interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] - out[j, i, jj] = out[j, i, jj] + - interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] + # Check if neighbor has corresponding child and if yes, establish connectivity + adjacent = adjacent_child(child, direction) + if has_child(t, neighbor_id, adjacent) + neighbor_child_id = t.child_ids[adjacent, neighbor_id] + opposite = opposite_direction(direction) + + t.neighbor_ids[direction, child_id] = neighbor_child_id + t.neighbor_ids[opposite, neighbor_child_id] = child_id end end - # apply inv(M) - @turbo for k in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d), - i in Base.OneTo(nnodes_1d) + return nothing +end + +# Refine given cells without rebalancing tree. +# +# Note: After a call to this method the tree may be unbalanced! +function refine_unbalanced!(t::AbstractTree, cell_ids, + sorted_unique_cell_ids = sort(unique(cell_ids))) + # Store actual ids refined cells (shifted due to previous insertions) + refined = zeros(Int, length(cell_ids)) + + # Loop over all cells that are to be refined + for (count, original_cell_id) in enumerate(sorted_unique_cell_ids) + # Determine actual cell id, taking into account previously inserted cells + n_children = n_children_per_cell(t) + cell_id = original_cell_id + (count - 1) * n_children + refined[count] = cell_id + + @assert !has_children(t, cell_id) "Non-leaf cell $cell_id cannot be refined" + + # Insert new cells directly behind parent (depth-first) + insert!(t, cell_id + 1, n_children) + + # Flip sign of refined cell such that we can easily find it later + t.original_cell_ids[cell_id] = -t.original_cell_ids[cell_id] + + # Initialize child cells (except neighbors) + for child in 1:n_children + child_id = cell_id + child + init_child!(t, cell_id, child, child_id) + end - out[i, j, k] = out[i, j, k] * inv_volume_weights_1d[i] * - inv_volume_weights_1d[j] * inv_volume_weights_1d[k] + # Initialize child cells (only neighbors) + # This separate loop is required since init_child_neighbors requires initialized parent-child + # relationships + for child in 1:n_children + child_id = cell_id + child + init_child_neighbors!(t, cell_id, child, child_id) + end end + + return refined end -# For now, this is mostly the same as `create_cache` for DGMultiFluxDiff{<:Polynomial}. -# In the future, we may modify it so that we can specialize additional parts of GaussSBP() solvers. -function create_cache(mesh::DGMultiMesh, equations, - dg::DGMultiFluxDiff{<:GaussSBP, <:Union{Quad, Hex}}, RealT, - uEltype) - - # call general Polynomial flux differencing constructor - cache = invoke(create_cache, - Tuple{typeof(mesh), typeof(equations), - DGMultiFluxDiff, typeof(RealT), typeof(uEltype)}, - mesh, equations, dg, RealT, uEltype) - - rd = dg.basis - @unpack md = mesh - - # for change of basis prior to the volume integral and entropy projection - r1D, _ = StartUpDG.gauss_lobatto_quad(0, 0, polydeg(dg)) - rq1D, _ = StartUpDG.gauss_quad(0, 0, polydeg(dg)) - interp_matrix_lobatto_to_gauss_1D = polynomial_interpolation_matrix(r1D, rq1D) - interp_matrix_gauss_to_lobatto_1D = polynomial_interpolation_matrix(rq1D, r1D) - NDIMS = ndims(rd.element_type) - interp_matrix_lobatto_to_gauss = SimpleKronecker(NDIMS, - interp_matrix_lobatto_to_gauss_1D, - uEltype) - interp_matrix_gauss_to_lobatto = SimpleKronecker(NDIMS, - interp_matrix_gauss_to_lobatto_1D, - uEltype) - inv_gauss_weights = inv.(rd.wq) - - # specialized operators to perform tensor product interpolation to faces for Gauss nodes - interp_matrix_gauss_to_face = TensorProductGaussFaceOperator(Interpolation(), dg) - projection_matrix_gauss_to_face = TensorProductGaussFaceOperator(Projection{Static.False()}(), - dg) - - # `LIFT` matrix for Gauss nodes - this is equivalent to `projection_matrix_gauss_to_face` scaled by `diagm(rd.wf)`, - # where `rd.wf` are Gauss node face quadrature weights. - gauss_LIFT = TensorProductGaussFaceOperator(Projection{Static.True()}(), dg) - - nvars = nvariables(equations) - rhs_volume_local_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) - for _ in 1:Threads.nthreads()] - gauss_volume_local_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) - for _ in 1:Threads.nthreads()] - - return (; cache..., projection_matrix_gauss_to_face, gauss_LIFT, inv_gauss_weights, - rhs_volume_local_threaded, gauss_volume_local_threaded, - interp_matrix_lobatto_to_gauss, interp_matrix_gauss_to_lobatto, - interp_matrix_gauss_to_face, - create_cache(mesh, equations, dg.volume_integral, dg, RealT, uEltype)...) # add cache specialized on the volume integral +# Refine entire tree by one level +function refine!(t::AbstractTree) + cells = @trixi_timeit timer() "collect all leaf cells" leaf_cells(t) + @trixi_timeit timer() "refine!" refine!(t, cells, cells) end -# by default, return an empty tuple for volume integral caches -create_cache(mesh, equations, volume_integral, dg, RealT, uEltype) = NamedTuple() - -# TODO: DGMulti. Address hard-coding of `entropy2cons!` and `cons2entropy!` for this function. -function entropy_projection!(cache, u, mesh::DGMultiMesh, equations, - dg::DGMultiFluxDiff{<:GaussSBP}) - rd = dg.basis - @unpack Vq = rd - @unpack VhP, entropy_var_values, u_values = cache - @unpack projected_entropy_var_values, entropy_projected_u_values = cache - @unpack interp_matrix_lobatto_to_gauss, interp_matrix_gauss_to_face = cache - - @threaded for e in eachelement(mesh, dg, cache) - apply_to_each_field(mul_by!(interp_matrix_lobatto_to_gauss), - view(u_values, :, e), view(u, :, e)) +# Refine given cells and rebalance tree. +# +# Note 1: Rebalancing is iterative, i.e., neighboring cells are refined if +# otherwise the 2:1 rule would be violated, which can cause more +# refinements. +# Note 2: Rebalancing currently only considers *Cartesian* neighbors, not diagonal neighbors! +function refine!(t::AbstractTree, cell_ids, + sorted_unique_cell_ids = sort(unique(cell_ids))) + # Reset original cell ids such that each cell knows its current id + reset_original_cell_ids!(t) + + # Refine all requested cells + refined = @trixi_timeit timer() "refine_unbalanced!" refine_unbalanced!(t, cell_ids, + sorted_unique_cell_ids) + refinement_count = length(refined) + + # Iteratively rebalance the tree until it does not change anymore + while length(refined) > 0 + refined = @trixi_timeit timer() "rebalance!" rebalance!(t, refined) + refinement_count += length(refined) end - # transform quadrature values to entropy variables - cons2entropy!(entropy_var_values, u_values, equations) + # Determine list of *original* cell ids that were refined + # Note: original_cell_ids contains the cell_id *before* refinement. At + # refinement, the refined cell's original_cell_ids value has its sign flipped + # to easily find it now. + refined_original_cells = @views(-t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0]) - volume_indices = Base.OneTo(rd.Nq) - face_indices = (rd.Nq + 1):(rd.Nq + rd.Nfq) + # Check if count of refinement cells matches information in original_cell_ids + @assert refinement_count==length(refined_original_cells) ("Mismatch in number of refined cells") - # Interpolate volume Gauss nodes to Gauss face nodes (note the layout of - # `projected_entropy_var_values = [vol pts; face pts]`). - entropy_var_face_values = view(projected_entropy_var_values, face_indices, :) - apply_to_each_field(mul_by!(interp_matrix_gauss_to_face), entropy_var_face_values, - entropy_var_values) + return refined_original_cells +end - # directly copy over volume values (no entropy projection required) - entropy_projected_volume_values = view(entropy_projected_u_values, volume_indices, - :) - @threaded for i in eachindex(u_values) - entropy_projected_volume_values[i] = u_values[i] +# Refine all leaf cells with coordinates in a given rectangular box +function refine_box!(t::AbstractTree{NDIMS}, coordinates_min, + coordinates_max) where {NDIMS} + for dim in 1:NDIMS + @assert coordinates_min[dim] cell_coordinates(t, cell_id))) + end - return nothing + # Refine cells + refine!(t, cells) +end + +# Convenience method for 1D +function refine_box!(t::AbstractTree{1}, coordinates_min::Real, coordinates_max::Real) + return refine_box!(t, [convert(Float64, coordinates_min)], + [convert(Float64, coordinates_max)]) end -# Assumes cache.flux_face_values is already computed. -# Enables tensor product evaluation of `LIFT isa TensorProductGaussFaceOperator`. -function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations, - surface_integral::SurfaceIntegralWeakForm, - dg::DGMultiFluxDiff{<:GaussSBP}, cache) - (; gauss_LIFT, gauss_volume_local_threaded) = cache +# Refine all leaf cells with coordinates in a given sphere +function refine_sphere!(t::AbstractTree{NDIMS}, center::SVector{NDIMS}, + radius) where {NDIMS} + @assert radius>=0 "Radius must be positive." - @threaded for e in eachelement(mesh, dg, cache) + # Find all leaf cells within sphere + cells = filter_leaf_cells(t) do cell_id + return sum(abs2, cell_coordinates(t, cell_id) - center) < radius^2 + end - # applies LIFT matrix, output is stored at Gauss nodes - gauss_volume_local = gauss_volume_local_threaded[Threads.threadid()] - apply_to_each_field(mul_by!(gauss_LIFT), gauss_volume_local, - view(cache.flux_face_values, :, e)) + # Refine cells + refine!(t, cells) +end - for i in eachindex(gauss_volume_local) - du[i, e] = du[i, e] + gauss_volume_local[i] +# Convenience function to allow passing center as a tuple +function refine_sphere!(t::AbstractTree{NDIMS}, center::NTuple{NDIMS}, + radius) where {NDIMS} + refine_sphere!(t, SVector(center), radius) +end + +# For the given cell ids, check if neighbors need to be refined to restore a rebalanced tree. +# +# Note 1: Rebalancing currently only considers *Cartesian* neighbors, not diagonal neighbors! +# Note 2: The current algorithm assumes that a previous refinement step has +# created level differences of at most 2. That is, before the previous +# refinement step, the tree was balanced. +function rebalance!(t::AbstractTree, refined_cell_ids) + # Create buffer for newly refined cells + to_refine = zeros(Int, n_directions(t) * length(refined_cell_ids)) + count = 0 + + # Iterate over cell ids that have previously been refined + for cell_id in refined_cell_ids + # Go over all potential neighbors of child cell + for direction in eachdirection(t) + # Continue if refined cell has a neighbor in that direction + if has_neighbor(t, cell_id, direction) + continue + end + + # Continue if refined cell has no coarse neighbor, since that would + # mean it there is no neighbor in that direction at all (domain + # boundary) + if !has_coarse_neighbor(t, cell_id, direction) + continue + end + + # Otherwise, the coarse neighbor exists and is not refined, thus it must + # be marked for refinement + coarse_neighbor_id = t.neighbor_ids[direction, t.parent_ids[cell_id]] + count += 1 + to_refine[count] = coarse_neighbor_id end end + + # Finally, refine all marked cells... + refined = refine_unbalanced!(t, unique(to_refine[1:count])) + + # ...and return list of refined cells + return refined end -@inline function flux_differencing_kernel!(du, u, element, mesh::DGMultiMesh, - have_nonconservative_terms, equations, - volume_flux, dg::DGMultiFluxDiff{<:GaussSBP}, - cache, alpha = true) - fluxdiff_local = cache.fluxdiff_local_threaded[Threads.threadid()] - fill!(fluxdiff_local, zero(eltype(fluxdiff_local))) - u_local = view(cache.entropy_projected_u_values, :, element) - - local_flux_differencing!(fluxdiff_local, u_local, element, - have_nonconservative_terms, - volume_flux, has_sparse_operators(dg), - mesh, equations, dg, cache) - - # convert `fluxdiff_local::Vector{<:SVector}` to `rhs_local::StructArray{<:SVector}` - # for faster performance when using `apply_to_each_field`. - rhs_local = cache.rhs_local_threaded[Threads.threadid()] - for i in Base.OneTo(length(fluxdiff_local)) - rhs_local[i] = fluxdiff_local[i] +# Refine given cells without rebalancing tree. +# +# Note: After a call to this method the tree may be unbalanced! +# function refine_unbalanced!(t::AbstractTree, cell_ids) end + +# Wrap single-cell refinements such that `sort(...)` does not complain +refine_unbalanced!(t::AbstractTree, cell_id::Int) = refine_unbalanced!(t, [cell_id]) + +# Coarsen entire tree by one level +function coarsen!(t::AbstractTree) + # Special case: if there is only one cell (root), there is nothing to do + if length(t) == 1 + return Int[] end - project_rhs_to_gauss_nodes!(du, rhs_local, element, mesh, dg, cache, alpha) + # Get list of unique parent ids for all leaf cells + parent_ids = unique(t.parent_ids[leaf_cells(t)]) + coarsen!(t, parent_ids) end -function project_rhs_to_gauss_nodes!(du, rhs_local, element, mesh::DGMultiMesh, - dg::DGMulti, cache, alpha = true) - - # Here, we exploit that under a Gauss nodal basis the structure of the projection - # matrix `Ph = [diagm(1 ./ wq), projection_matrix_gauss_to_face]` such that - # `Ph * [u; uf] = (u ./ wq) + projection_matrix_gauss_to_face * uf`. - volume_indices = Base.OneTo(dg.basis.Nq) - face_indices = (dg.basis.Nq + 1):(dg.basis.Nq + dg.basis.Nfq) - local_volume_flux = view(rhs_local, volume_indices) - local_face_flux = view(rhs_local, face_indices) - - # initialize rhs_volume_local = projection_matrix_gauss_to_face * local_face_flux - rhs_volume_local = cache.rhs_volume_local_threaded[Threads.threadid()] - apply_to_each_field(mul_by!(cache.projection_matrix_gauss_to_face), - rhs_volume_local, local_face_flux) - - # accumulate volume contributions at Gauss nodes - for i in eachindex(rhs_volume_local) - du_local = rhs_volume_local[i] + - local_volume_flux[i] * cache.inv_gauss_weights[i] - du[i, element] = du[i, element] + alpha * du_local +# Coarsen given *parent* cells (= these cells must have children who are all +# leaf cells) while retaining a balanced tree. +# +# A cell to be coarsened might cause an unbalanced tree if the neighboring cell +# was already refined. Since it is generally not desired that cells are +# coarsened without specifically asking for it, these cells will then *not* be +# coarsened. +function coarsen!(t::AbstractTree, cell_ids::AbstractArray{Int}) + # Return early if array is empty + if length(cell_ids) == 0 + return Int[] end -end -function calc_volume_integral!(du, u, mesh::DGMultiMesh, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralFluxDifferencing, - dg::DGMultiFluxDiff{<:GaussSBP}, cache) - @threaded for e in eachelement(mesh, dg, cache) - flux_differencing_kernel!(du, u, e, mesh, - have_nonconservative_terms, equations, - volume_integral.volume_flux, dg, cache) + # Reset original cell ids such that each cell knows its current id + reset_original_cell_ids!(t) + + # To maximize the number of cells that may be coarsened, start with the cells at the highest level + sorted_by_level = sort(cell_ids, by = i -> t.levels[i]) + + # Keep track of number of cells that were actually coarsened + n_coarsened = 0 + + # Local function to adjust cell ids after some cells have been removed + function adjust_cell_ids!(cell_ids, coarsened_cell_id, count) + for (id, cell_id) in enumerate(cell_ids) + if cell_id > coarsened_cell_id + cell_ids[id] = cell_id - count + end + end end -end -# interpolate back to Lobatto nodes after applying the inverse Jacobian at Gauss points -function invert_jacobian_and_interpolate!(du, mesh::DGMultiMesh, equations, - dg::DGMultiFluxDiff{<:GaussSBP}, cache; - scaling = -1) - (; interp_matrix_gauss_to_lobatto, rhs_volume_local_threaded, invJ) = cache + # Iterate backwards over cells to coarsen + while true + # Retrieve next cell or quit + if length(sorted_by_level) > 0 + coarse_cell_id = pop!(sorted_by_level) + else + break + end - @threaded for e in eachelement(mesh, dg, cache) - rhs_volume_local = rhs_volume_local_threaded[Threads.threadid()] + # Ensure that cell has children (violation is an error) + if !has_children(t, coarse_cell_id) + error("cell is leaf and cannot be coarsened to: $coarse_cell_id") + end + + # Ensure that all child cells are leaf cells (violation is an error) + for child in 1:n_children_per_cell(t) + if has_child(t, coarse_cell_id, child) + if !is_leaf(t, t.child_ids[child, coarse_cell_id]) + error("cell $coarse_cell_id has child cell at position $child that is not a leaf cell") + end + end + end - # At this point, `rhs_volume_local` should still be stored at Gauss points. - # We scale it by the inverse Jacobian before transforming back to Lobatto. - for i in eachindex(rhs_volume_local) - rhs_volume_local[i] = du[i, e] * invJ[i, e] * scaling + # Check if coarse cell has refined neighbors that would prevent coarsening + skip = false + # Iterate over all children (which are to be removed) + for child in 1:n_children_per_cell(t) + # Continue if child does not exist + if !has_child(t, coarse_cell_id, child) + continue + end + child_id = t.child_ids[child, coarse_cell_id] + + # Go over all neighbors of child cell. If it has a neighbor that is *not* + # a sibling and that is not a leaf cell, we cannot coarsen its parent + # without creating an unbalanced tree. + for direction in eachdirection(t) + # Continue if neighbor would be a sibling + if has_sibling(child, direction) + continue + end + + # Continue if child cell has no neighbor in that direction + if !has_neighbor(t, child_id, direction) + continue + end + neighbor_id = t.neighbor_ids[direction, child_id] + + if !has_children(t, neighbor_id) + continue + end + + # If neighbor is not a sibling, is existing, and has children, do not coarsen + skip = true + break + end + end + # Skip if a neighboring cell prevents coarsening + if skip + continue end - # Interpolate result back to Lobatto nodes for ease of analysis, visualization - apply_to_each_field(mul_by!(interp_matrix_gauss_to_lobatto), - view(du, :, e), rhs_volume_local) + # Flip sign of cell to be coarsened to such that we can easily find it + t.original_cell_ids[coarse_cell_id] = -t.original_cell_ids[coarse_cell_id] + + # If a coarse cell has children that are all leaf cells, they must follow + # immediately due to depth-first ordering of the tree + count = n_children(t, coarse_cell_id) + @assert count==n_children_per_cell(t) "cell $coarse_cell_id does not have all child cells" + remove_shift!(t, coarse_cell_id + 1, coarse_cell_id + count) + + # Take into account shifts in tree that alters cell ids + adjust_cell_ids!(sorted_by_level, coarse_cell_id, count) + + # Keep track of number of coarsened cells + n_coarsened += 1 end + + # Determine list of *original* cell ids that were coarsened to + # Note: original_cell_ids contains the cell_id *before* coarsening. At + # coarsening, the coarsened parent cell's original_cell_ids value has its sign flipped + # to easily find it now. + @views coarsened_original_cells = (-t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0]) + + # Check if count of coarsened cells matches information in original_cell_ids + @assert n_coarsened==length(coarsened_original_cells) ("Mismatch in number of coarsened cells") + + return coarsened_original_cells end -# Specialize RHS so that we can call `invert_jacobian_and_interpolate!` instead of just `invert_jacobian!`, -# since `invert_jacobian!` is also used in other places (e.g., parabolic terms). -function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions::BC, - source_terms::Source, dg::DGMultiFluxDiff{<:GaussSBP}, - cache) where {Source, BC} - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) - - # this function evaluates the solution at volume and face quadrature points (which was previously - # done in `prolong2interfaces` and `calc_volume_integral`) - @trixi_timeit timer() "entropy_projection!" begin - entropy_projection!(cache, u, mesh, equations, dg) - end +# Wrap single-cell coarsening such that `sort(...)` does not complain +coarsen!(t::AbstractTree, cell_id::Int) = coarsen!(t::AbstractTree, [cell_id]) - # `du` is stored at Gauss nodes here - @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, u, mesh, - have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) +# Coarsen all viable parent cells with coordinates in a given rectangular box +function coarsen_box!(t::AbstractTree{NDIMS}, coordinates_min::AbstractArray{Float64}, + coordinates_max::AbstractArray{Float64}) where {NDIMS} + for dim in 1:NDIMS + @assert coordinates_min[dim] cell_coordinates(t, cell_id))) end - @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux!(cache, t, boundary_conditions, mesh, - have_nonconservative_terms(equations), equations, dg) - end + # Get list of unique parent ids for all leaf cells + parent_ids = unique(t.parent_ids[leaves]) - # `du` is stored at Gauss nodes here - @trixi_timeit timer() "surface integral" begin - calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache) + # Filter parent ids to be within box + parents = filter(parent_ids) do cell_id + return (all(coordinates_min .< cell_coordinates(t, cell_id)) && + all(coordinates_max .> cell_coordinates(t, cell_id))) end - # invert Jacobian and map `du` from Gauss to Lobatto nodes - @trixi_timeit timer() "Jacobian" begin - invert_jacobian_and_interpolate!(du, mesh, equations, dg, cache) - end + # Coarsen cells + coarsen!(t, parents) +end - @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, mesh, equations, dg, cache) +# Convenience method for 1D +function coarsen_box!(t::AbstractTree{1}, coordinates_min::Real, coordinates_max::Real) + return coarsen_box!(t, [convert(Float64, coordinates_min)], + [convert(Float64, coordinates_max)]) +end + +# Return coordinates of a child cell based on its relative position to the parent. +function child_coordinates(::AbstractTree{NDIMS}, parent_coordinates, + parent_length::Number, child::Int) where {NDIMS} + # Calculate length of child cells + child_length = parent_length / 2 + return SVector(ntuple(d -> parent_coordinates[d] + + child_sign(child, d) * child_length / 2, Val(NDIMS))) +end + +# Reset range of cells to values that are prone to cause errors as soon as they are used. +# +# Rationale: If an invalid cell is accidentally used, we want to know it as soon as possible. +# function invalidate!(t::AbstractTree, first::Int, last::Int) end +invalidate!(t::AbstractTree, id::Int) = invalidate!(t, id, id) +invalidate!(t::AbstractTree) = invalidate!(t, 1, length(t)) + +# Delete connectivity with parents/children/neighbors before cells are erased +function delete_connectivity!(t::AbstractTree, first::Int, last::Int) + @assert first > 0 + @assert first <= last + @assert last <= t.capacity + 1 + + # Iterate over all cells + for cell_id in first:last + # Delete connectivity from parent cell + if has_parent(t, cell_id) + parent_id = t.parent_ids[cell_id] + for child in 1:n_children_per_cell(t) + if t.child_ids[child, parent_id] == cell_id + t.child_ids[child, parent_id] = 0 + break + end + end + end + + # Delete connectivity from child cells + for child in 1:n_children_per_cell(t) + if has_child(t, cell_id, child) + t.parent_ids[t._child_ids[child, cell_id]] = 0 + end + end + + # Delete connectivity from neighboring cells + for direction in eachdirection(t) + if has_neighbor(t, cell_id, direction) + t.neighbor_ids[opposite_direction(direction), t.neighbor_ids[direction, cell_id]] = 0 + end + end end +end - return nothing +# Move connectivity with parents/children/neighbors after cells have been moved +function move_connectivity!(t::AbstractTree, first::Int, last::Int, destination::Int) + @assert first > 0 + @assert first <= last + @assert last <= t.capacity + 1 + @assert destination > 0 + @assert destination <= t.capacity + 1 + + # Strategy + # 1) Loop over moved cells (at target location) + # 2) Check if parent/children/neighbors connections are to a cell that was moved + # a) if cell was moved: apply offset to current cell + # b) if cell was not moved: go to connected cell and update connectivity there + + offset = destination - first + has_moved(n) = (first <= n <= last) + + for source in first:last + target = source + offset + + # Update parent + if has_parent(t, target) + # Get parent cell + parent_id = t.parent_ids[target] + if has_moved(parent_id) + # If parent itself was moved, just update parent id accordingly + t.parent_ids[target] += offset + else + # If parent was not moved, update its corresponding child id + for child in 1:n_children_per_cell(t) + if t.child_ids[child, parent_id] == source + t.child_ids[child, parent_id] = target + end + end + end + end + + # Update children + for child in 1:n_children_per_cell(t) + if has_child(t, target, child) + # Get child cell + child_id = t.child_ids[child, target] + if has_moved(child_id) + # If child itself was moved, just update child id accordingly + t.child_ids[child, target] += offset + else + # If child was not moved, update its parent id + t.parent_ids[child_id] = target + end + end + end + + # Update neighbors + for direction in eachdirection(t) + if has_neighbor(t, target, direction) + # Get neighbor cell + neighbor_id = t.neighbor_ids[direction, target] + if has_moved(neighbor_id) + # If neighbor itself was moved, just update neighbor id accordingly + t.neighbor_ids[direction, target] += offset + else + # If neighbor was not moved, update its opposing neighbor id + t.neighbor_ids[opposite_direction(direction), neighbor_id] = target + end + end + end + end end + +# Raw copy operation for ranges of cells. +# +# This method is used by the higher-level copy operations for AbstractContainer +# function raw_copy!(target::AbstractTree, source::AbstractTree, first::Int, last::Int, destination::Int) end + +# Reset data structures by recreating all internal storage containers and invalidating all elements +# function reset_data_structures!(t::AbstractTree{NDIMS}) where NDIMS end + end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index 43649a5dc6e..eaa6ab5cee2 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -49,7 +49,7 @@ function start_mpi_send!(mpi_cache::P4estMPICache, mesh, equations, dg, cache) data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1) n_small_elements = 2^(ndims(mesh) - 1) - for d in eachindex(mpi_cache.mpi_neighbor_ranks) + for d in 1:length(mpi_cache.mpi_neighbor_ranks) send_buffer = mpi_cache.mpi_send_buffers[d] for (index, interface) in enumerate(mpi_cache.mpi_neighbor_interfaces[d]) @@ -505,7 +505,7 @@ function exchange_normal_directions!(mpi_mortars, mpi_cache, # Create buffers and requests send_buffers = Vector{Vector{RealT}}(undef, length(mpi_neighbor_mortars)) recv_buffers = Vector{Vector{RealT}}(undef, length(mpi_neighbor_mortars)) - for index in eachindex(mpi_neighbor_mortars) + for index in 1:length(mpi_neighbor_mortars) send_buffers[index] = Vector{RealT}(undef, length(mpi_neighbor_mortars[index]) * n_small_elements * data_size) @@ -519,7 +519,7 @@ function exchange_normal_directions!(mpi_mortars, mpi_cache, recv_requests = Vector{MPI.Request}(undef, length(mpi_neighbor_mortars)) # Fill send buffers - for d in eachindex(mpi_neighbor_ranks) + for d in 1:length(mpi_neighbor_ranks) send_buffer = send_buffers[d] for (index, mortar) in enumerate(mpi_neighbor_mortars[d]) diff --git a/src/solvers/dgsem_tree/dg_parallel.jl b/src/solvers/dgsem_tree/dg_parallel.jl index f1e78b2b64a..c614fe0d0e6 100644 --- a/src/solvers/dgsem_tree/dg_parallel.jl +++ b/src/solvers/dgsem_tree/dg_parallel.jl @@ -13,7 +13,7 @@ function init_mpi_data_structures(mpi_neighbor_interfaces, mpi_neighbor_mortars, n_small_elements = 2^(n_dims - 1) mpi_send_buffers = Vector{Vector{uEltype}}(undef, length(mpi_neighbor_interfaces)) mpi_recv_buffers = Vector{Vector{uEltype}}(undef, length(mpi_neighbor_interfaces)) - for index in eachindex(mpi_neighbor_interfaces) + for index in 1:length(mpi_neighbor_interfaces) mpi_send_buffers[index] = Vector{uEltype}(undef, length(mpi_neighbor_interfaces[index]) * data_size + From 2e79eaa663cb4dd3e39f7badb1400c1d0a3ad13d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 8 Apr 2024 09:31:34 +0200 Subject: [PATCH 4/9] revert --- .../dgmulti/flux_differencing_gauss_sbp.jl | 1162 +++++++---------- 1 file changed, 488 insertions(+), 674 deletions(-) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index 469189ff50c..7c9227b5c3c 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -1,779 +1,593 @@ -# 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 -abstract type AbstractTree{NDIMS} <: AbstractContainer end +# ========= GaussSBP approximation types ============ +# Note: we define type aliases outside of the @muladd block to avoid Revise breaking when code +# inside the @muladd block is edited. See https://github.com/trixi-framework/Trixi.jl/issues/801 +# for more details. -# Type traits to obtain dimension -@inline Base.ndims(::AbstractTree{NDIMS}) where {NDIMS} = NDIMS +# `GaussSBP` is a type alias for a StartUpDG type (e.g., Gauss nodes on quads/hexes) +const GaussSBP = Polynomial{Gauss} -# Auxiliary methods to allow semantic queries on the tree -# Check whether cell has parent cell -has_parent(t::AbstractTree, cell_id::Int) = t.parent_ids[cell_id] > 0 - -# Count number of children for a given cell -function n_children(t::AbstractTree, cell_id::Int) - count(x -> (x > 0), @view t.child_ids[:, cell_id]) +function tensor_product_quadrature(element_type::Line, r1D, w1D) + return r1D, w1D end -# Check whether cell has any child cell -has_children(t::AbstractTree, cell_id::Int) = n_children(t, cell_id) > 0 - -# Check whether cell is leaf cell -is_leaf(t::AbstractTree, cell_id::Int) = !has_children(t, cell_id) - -# Check whether cell has specific child cell -has_child(t::AbstractTree, cell_id::Int, child::Int) = t.child_ids[child, cell_id] > 0 - -# Check if cell has a neighbor at the same refinement level in the given direction -function has_neighbor(t::AbstractTree, cell_id::Int, direction::Int) - t.neighbor_ids[direction, cell_id] > 0 -end - -# Check if cell has a coarse neighbor, i.e., with one refinement level lower -function has_coarse_neighbor(t::AbstractTree, cell_id::Int, direction::Int) - return has_parent(t, cell_id) && has_neighbor(t, t.parent_ids[cell_id], direction) +function tensor_product_quadrature(element_type::Quad, r1D, w1D) + sq, rq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D)) + ws, wr = vec.(StartUpDG.NodesAndModes.meshgrid(w1D)) + wq = wr .* ws + return rq, sq, wq end -# Check if cell has any neighbor (same-level or lower-level) -function has_any_neighbor(t::AbstractTree, cell_id::Int, direction::Int) - return has_neighbor(t, cell_id, direction) || - has_coarse_neighbor(t, cell_id, direction) +function tensor_product_quadrature(element_type::Hex, r1D, w1D) + rq, sq, tq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D, r1D, r1D)) + wr, ws, wt = vec.(StartUpDG.NodesAndModes.meshgrid(w1D, w1D, w1D)) + wq = wr .* ws .* wt + return rq, sq, tq, wq end -# Check if cell is own cell, i.e., belongs to this MPI rank -is_own_cell(t::AbstractTree, cell_id) = true - -# Return cell length for a given level -length_at_level(t::AbstractTree, level::Int) = t.length_level_0 / 2^level - -# Return cell length for a given cell -length_at_cell(t::AbstractTree, cell_id::Int) = length_at_level(t, t.levels[cell_id]) - -# Return minimum level of any leaf cell -minimum_level(t::AbstractTree) = minimum(t.levels[leaf_cells(t)]) +# type parameters for `TensorProductFaceOperator`. +abstract type AbstractGaussOperator end +struct Interpolation <: AbstractGaussOperator end +# - `Projection{ScaleByFaceWeights=Static.False()}` corresponds to the operator `projection_matrix_gauss_to_face = M \ Vf'`, +# which is used in `VolumeIntegralFluxDifferencing`. +# - `Projection{ScaleByFaceWeights=Static.True()}` corresponds to the quadrature-based lifting +# operator `LIFT = M \ (Vf' * diagm(rd.wf))`, which is used in `SurfaceIntegralWeakForm` +struct Projection{ScaleByFaceWeights} <: AbstractGaussOperator end -# Return maximum level of any leaf cell -maximum_level(t::AbstractTree) = maximum(t.levels[leaf_cells(t)]) +# used to dispatch for different Gauss interpolation operators +abstract type AbstractTensorProductGaussOperator end -# Check if tree is periodic -isperiodic(t::AbstractTree) = all(t.periodicity) -isperiodic(t::AbstractTree, dimension) = t.periodicity[dimension] - -# Auxiliary methods for often-required calculations -# Number of potential child cells -n_children_per_cell(::AbstractTree{NDIMS}) where {NDIMS} = 2^NDIMS - -# Number of directions -# -# Directions are indicated by numbers from 1 to 2*ndims: -# 1 -> -x -# 2 -> +x -# 3 -> -y -# 4 -> +y -# 5 -> -z -# 6 -> +z -@inline n_directions(::AbstractTree{NDIMS}) where {NDIMS} = 2 * NDIMS -# TODO: Taal performance, 1:n_directions(tree) vs. Base.OneTo(n_directions(tree)) vs. SOneTo(n_directions(tree)) -""" - eachdirection(tree::AbstractTree) - -Return an iterator over the indices that specify the location in relevant data structures -for the directions in `AbstractTree`. -In particular, not the directions themselves are returned. -""" -@inline eachdirection(tree::AbstractTree) = Base.OneTo(n_directions(tree)) - -# For a given direction, return its opposite direction -# -# dir -> opp -# 1 -> 2 -# 2 -> 1 -# 3 -> 4 -# 4 -> 3 -# 5 -> 6 -# 6 -> 5 -opposite_direction(direction::Int) = direction + 1 - 2 * ((direction + 1) % 2) - -# For a given child position (from 1 to 8) and dimension (from 1 to 3), -# calculate a child cell's position relative to its parent cell. +# TensorProductGaussFaceOperator{Tmat, Ti} # -# Essentially calculates the following -# dim=1 dim=2 dim=3 -# child x y z -# 1 - - - -# 2 + - - -# 3 - + - -# 4 + + - -# 5 - - + -# 6 + - + -# 7 - + + -# 8 + + + -# child_sign(child::Int, dim::Int) = 1 - 2 * (div(child + 2^(dim - 1) - 1, 2^(dim-1)) % 2) -# Since we use only a fixed number of dimensions, we use a lookup table for improved performance. -const _child_signs = [-1 -1 -1; - +1 -1 -1; - -1 +1 -1; - +1 +1 -1; - -1 -1 +1; - +1 -1 +1; - -1 +1 +1; - +1 +1 +1] -child_sign(child::Int, dim::Int) = _child_signs[child, dim] - -# For each child position (1 to 8) and a given direction (from 1 to 6), return -# neighboring child position. -const _adjacent_child_ids = [2 2 3 3 5 5; - 1 1 4 4 6 6; - 4 4 1 1 7 7; - 3 3 2 2 8 8; - 6 6 7 7 1 1; - 5 5 8 8 2 2; - 8 8 5 5 3 3; - 7 7 6 6 4 4] -adjacent_child(child::Int, direction::Int) = _adjacent_child_ids[child, direction] - -# For each child position (1 to 8) and a given direction (from 1 to 6), return -# if neighbor is a sibling -function has_sibling(child::Int, direction::Int) - return (child_sign(child, div(direction + 1, 2)) * (-1)^(direction - 1)) > 0 +# Data for performing tensor product interpolation from volume nodes to face nodes. +struct TensorProductGaussFaceOperator{NDIMS, OperatorType <: AbstractGaussOperator, + Tmat, Tweights, Tfweights, Tindices} <: + AbstractTensorProductGaussOperator + interp_matrix_gauss_to_face_1d::Tmat + inv_volume_weights_1d::Tweights + face_weights::Tfweights + face_indices_tensor_product::Tindices + nnodes_1d::Int + nfaces::Int end -# Obtain leaf cells that fulfill a given criterion. -# -# The function `f` is passed the cell id of each leaf cell -# as an argument. -function filter_leaf_cells(f, t::AbstractTree) - filtered = Vector{Int}(undef, length(t)) - count = 0 - for cell_id in 1:length(t) - if is_leaf(t, cell_id) && f(cell_id) - count += 1 - filtered[count] = cell_id - end - end +# constructor for a 2D operator +function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, + dg::DGMulti{2, Quad, GaussSBP}) + rd = dg.basis - return filtered[1:count] -end + rq1D, wq1D = StartUpDG.gauss_quad(0, 0, polydeg(dg)) + interp_matrix_gauss_to_face_1d = polynomial_interpolation_matrix(rq1D, [-1; 1]) -# Return an array with the ids of all leaf cells -leaf_cells(t::AbstractTree) = filter_leaf_cells((cell_id) -> true, t) + nnodes_1d = length(rq1D) -# Return an array with the ids of all leaf cells for a given rank -leaf_cells_by_rank(t::AbstractTree, rank) = leaf_cells(t) - -# Return an array with the ids of all local leaf cells -local_leaf_cells(t::AbstractTree) = leaf_cells(t) - -# Count the number of leaf cells. -count_leaf_cells(t::AbstractTree) = length(leaf_cells(t)) + # Permutation of indices in a tensor product form + num_faces = StartUpDG.num_faces(rd.element_type) + indices = reshape(1:length(rd.rf), nnodes_1d, num_faces) + face_indices_tensor_product = zeros(Int, 2, nnodes_1d, ndims(rd.element_type)) + for i in 1:nnodes_1d # loop over nodes in one face + face_indices_tensor_product[:, i, 1] .= indices[i, 1:2] + face_indices_tensor_product[:, i, 2] .= indices[i, 3:4] + end -@inline function cell_coordinates(t::AbstractTree{NDIMS}, cell) where {NDIMS} - SVector(ntuple(d -> t.coordinates[d, cell], Val(NDIMS))) + T_op = typeof(operator) + Tm = typeof(interp_matrix_gauss_to_face_1d) + Tw = typeof(inv.(wq1D)) + Tf = typeof(rd.wf) + Ti = typeof(face_indices_tensor_product) + return TensorProductGaussFaceOperator{2, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d, + inv.(wq1D), rd.wf, + face_indices_tensor_product, + nnodes_1d, num_faces) end -@inline function set_cell_coordinates!(t::AbstractTree{NDIMS}, coords, - cell) where {NDIMS} - for d in 1:NDIMS - t.coordinates[d, cell] = coords[d] +# constructor for a 3D operator +function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, + dg::DGMulti{3, Hex, GaussSBP}) + rd = dg.basis + + rq1D, wq1D = StartUpDG.gauss_quad(0, 0, polydeg(dg)) + interp_matrix_gauss_to_face_1d = polynomial_interpolation_matrix(rq1D, [-1; 1]) + + nnodes_1d = length(rq1D) + + # Permutation of indices in a tensor product form + num_faces = StartUpDG.num_faces(rd.element_type) + indices = reshape(1:length(rd.rf), nnodes_1d, nnodes_1d, num_faces) + face_indices_tensor_product = zeros(Int, 2, nnodes_1d, nnodes_1d, + ndims(rd.element_type)) + for j in 1:nnodes_1d, i in 1:nnodes_1d # loop over nodes in one face + face_indices_tensor_product[:, i, j, 1] .= indices[i, j, 1:2] + face_indices_tensor_product[:, i, j, 2] .= indices[i, j, 3:4] + face_indices_tensor_product[:, i, j, 3] .= indices[i, j, 5:6] end -end -# Determine if point is located inside cell -function is_point_in_cell(t::AbstractTree, point_coordinates, cell_id) - cell_length = length_at_cell(t, cell_id) - cell_coordinates_ = cell_coordinates(t, cell_id) - min_coordinates = cell_coordinates_ .- cell_length / 2 - max_coordinates = cell_coordinates_ .+ cell_length / 2 + T_op = typeof(operator) + Tm = typeof(interp_matrix_gauss_to_face_1d) + Tw = typeof(inv.(wq1D)) + Tf = typeof(rd.wf) + Ti = typeof(face_indices_tensor_product) + return TensorProductGaussFaceOperator{3, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d, + inv.(wq1D), rd.wf, + face_indices_tensor_product, + nnodes_1d, num_faces) +end - return all(min_coordinates .<= point_coordinates .<= max_coordinates) +# specialize behavior of `mul_by!(A)` where `A isa TensorProductGaussFaceOperator)` +@inline function mul_by!(A::AbstractTensorProductGaussOperator) + return (out, x) -> tensor_product_gauss_face_operator!(out, A, x) end -# Store cell id in each cell to use for post-AMR analysis -function reset_original_cell_ids!(t::AbstractTree) - t.original_cell_ids[1:length(t)] .= 1:length(t) +@inline function tensor_product_gauss_face_operator!(out::AbstractMatrix, + A::AbstractTensorProductGaussOperator, + x::AbstractMatrix) + @threaded for col in Base.OneTo(size(out, 2)) + tensor_product_gauss_face_operator!(view(out, :, col), A, view(x, :, col)) + end end -# Efficiently perform uniform refinement up to a given level (works only on mesh with a single cell) -function refine_uniformly!(t::AbstractTree, max_level) - @assert length(t)==1 "efficient uniform refinement only works for a newly created tree" - @assert max_level>=0 "the uniform refinement level must be non-zero" +# 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 - # Calculate size of final tree and resize tree - total_length = 1 - for level in 1:max_level - total_length += n_children_per_cell(t)^level +#! format: off +# Interpolates values from volume Gauss nodes to face nodes on one element. +@inline function tensor_product_gauss_face_operator!(out::AbstractVector, + A::TensorProductGaussFaceOperator{2, Interpolation}, + x_in::AbstractVector) +#! format: on + (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A + (; nnodes_1d) = A + + fill!(out, zero(eltype(out))) + + # for 2D GaussSBP nodes, the indexing is first in x, then in y + x = reshape(x_in, nnodes_1d, nnodes_1d) + + # interpolation in the x-direction + @turbo for i in Base.OneTo(nnodes_1d) # loop over nodes in a face + index_left = face_indices_tensor_product[1, i, 1] + index_right = face_indices_tensor_product[2, i, 1] + for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes + out[index_left] = out[index_left] + + interp_matrix_gauss_to_face_1d[1, jj] * x[jj, i] + out[index_right] = out[index_right] + + interp_matrix_gauss_to_face_1d[2, jj] * x[jj, i] + end end - resize!(t, total_length) - # Traverse tree to set parent-child relationships - init_children!(t, 1, max_level) - - # Set all neighbor relationships - init_neighbors!(t, max_level) + # interpolation in the y-direction + @turbo for i in Base.OneTo(nnodes_1d) # loop over nodes in a face + index_left = face_indices_tensor_product[1, i, 2] + index_right = face_indices_tensor_product[2, i, 2] + for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes + out[index_left] = out[index_left] + + interp_matrix_gauss_to_face_1d[1, jj] * x[i, jj] + out[index_right] = out[index_right] + + interp_matrix_gauss_to_face_1d[2, jj] * x[i, jj] + end + end end -# Recursively initialize children up to level `max_level` in depth-first ordering, starting with -# cell `cell_id` and set all information except neighbor relations (see `init_neighbors!`). -# -# Return the number of offspring of the initialized cell plus one -function init_children!(t::AbstractTree, cell_id, max_level) - # Stop recursion if max_level has been reached - if t.levels[cell_id] >= max_level - return 1 - else - # Initialize each child cell, counting the total number of offspring - n_offspring = 0 - for child in 1:n_children_per_cell(t) - # Get cell id of child - child_id = cell_id + 1 + n_offspring - - # Initialize child cell (except neighbors) - init_child!(t, cell_id, child, child_id) - - # Recursively initialize child cell - n_offspring += init_children!(t, child_id, max_level) +# Interpolates values from volume Gauss nodes to face nodes on one element. +#! format: off +@inline function tensor_product_gauss_face_operator!(out::AbstractVector, + A::TensorProductGaussFaceOperator{3, Interpolation}, + x::AbstractVector) +#! format: on + (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A + (; nnodes_1d) = A + + fill!(out, zero(eltype(out))) + + # for 3D GaussSBP nodes, the indexing is first in y, then x, then z. + x = reshape(x, nnodes_1d, nnodes_1d, nnodes_1d) + + # interpolation in the y-direction + @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face + index_left = face_indices_tensor_product[1, i, j, 2] + index_right = face_indices_tensor_product[2, i, j, 2] + for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes + out[index_left] = out[index_left] + + interp_matrix_gauss_to_face_1d[1, jj] * x[jj, i, j] + out[index_right] = out[index_right] + + interp_matrix_gauss_to_face_1d[2, jj] * x[jj, i, j] end - - return n_offspring + 1 end -end -# Iteratively set all neighbor relations, starting at an initialized level 0 cell. Assume that -# parent-child relations have already been initialized (see `init_children!`). -function init_neighbors!(t::AbstractTree, max_level = maximum_level(t)) - @assert all(n >= 0 for n in t.neighbor_ids[:, 1]) "level 0 cell neighbors must be initialized" - - # Initialize neighbors level by level - for level in 1:max_level - # Walk entire tree, starting from level 0 cell - for cell_id in 1:length(t) - # Skip cells whose immediate children are already initialized *or* whose level is too high for this round - if t.levels[cell_id] != level - 1 - continue - end - - # Iterate over children and set neighbor information - for child in 1:n_children_per_cell(t) - child_id = t.child_ids[child, cell_id] - init_child_neighbors!(t, cell_id, child, child_id) - end + # interpolation in the x-direction + @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face + index_left = face_indices_tensor_product[1, i, j, 1] + index_right = face_indices_tensor_product[2, i, j, 1] + for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes + out[index_left] = out[index_left] + + interp_matrix_gauss_to_face_1d[1, jj] * x[i, jj, j] + out[index_right] = out[index_right] + + interp_matrix_gauss_to_face_1d[2, jj] * x[i, jj, j] end end - return nothing + # interpolation in the z-direction + @turbo for i in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d) # loop over nodes in a face + index_left = face_indices_tensor_product[1, i, j, 3] + index_right = face_indices_tensor_product[2, i, j, 3] + for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes + # The ordering (i,j) -> (j,i) needs to be reversed for this last face. + # This is due to way we define face nodes for Hex() types in StartUpDG.jl. + out[index_left] = out[index_left] + + interp_matrix_gauss_to_face_1d[1, jj] * x[j, i, jj] + out[index_right] = out[index_right] + + interp_matrix_gauss_to_face_1d[2, jj] * x[j, i, jj] + end + end end -# Initialize the neighbors of child cell `child_id` based on parent cell `cell_id` -function init_child_neighbors!(t::AbstractTree, cell_id, child, child_id) - t.neighbor_ids[:, child_id] .= zero(eltype(t.neighbor_ids)) - for direction in eachdirection(t) - # If neighbor is a sibling, establish one-sided connectivity - # Note: two-sided is not necessary, as each sibling will do this - if has_sibling(child, direction) - adjacent = adjacent_child(child, direction) - neighbor_id = t.child_ids[adjacent, cell_id] - - t.neighbor_ids[direction, child_id] = neighbor_id - continue +# Projects face node values to volume Gauss nodes on one element. +#! format: off +@inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector, + A::TensorProductGaussFaceOperator{2, Projection{ApplyFaceWeights}}, + x::AbstractVector) where {ApplyFaceWeights} +#! format: on + (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A + (; inv_volume_weights_1d, nnodes_1d) = A + + fill!(out_vec, zero(eltype(out_vec))) + + # As of Julia 1.9, Base.ReshapedArray does not produce allocations when setting values. + # Thus, Base.ReshapedArray should be used if you are setting values in the array. + # `reshape` is fine if you are only accessing values. + # Note that, for 2D GaussSBP nodes, the indexing is first in x, then y + out = Base.ReshapedArray(out_vec, (nnodes_1d, nnodes_1d), ()) + + if ApplyFaceWeights == true + @turbo for i in eachindex(x) + x[i] = x[i] * A.face_weights[i] end + end - # Skip if original cell does have no neighbor in direction - if !has_neighbor(t, cell_id, direction) - continue + # interpolation in the x-direction + @turbo for i in Base.OneTo(nnodes_1d) # loop over face nodes + index_left = face_indices_tensor_product[1, i, 1] + index_right = face_indices_tensor_product[2, i, 1] + for jj in Base.OneTo(nnodes_1d) # loop over a line of volume nodes + out[jj, i] = out[jj, i] + + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] + out[jj, i] = out[jj, i] + + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] end + end - # Otherwise, check if neighbor has children - if not, skip again - neighbor_id = t.neighbor_ids[direction, cell_id] - if !has_children(t, neighbor_id) - continue + # interpolation in the y-direction + @turbo for i in Base.OneTo(nnodes_1d) + index_left = face_indices_tensor_product[1, i, 2] + index_right = face_indices_tensor_product[2, i, 2] + # loop over a line of volume nodes + for jj in Base.OneTo(nnodes_1d) + out[i, jj] = out[i, jj] + + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] + out[i, jj] = out[i, jj] + + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] end + end - # Check if neighbor has corresponding child and if yes, establish connectivity - adjacent = adjacent_child(child, direction) - if has_child(t, neighbor_id, adjacent) - neighbor_child_id = t.child_ids[adjacent, neighbor_id] - opposite = opposite_direction(direction) + # apply inv(M) + @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) + out[i, j] = out[i, j] * inv_volume_weights_1d[i] * inv_volume_weights_1d[j] + end +end - t.neighbor_ids[direction, child_id] = neighbor_child_id - t.neighbor_ids[opposite, neighbor_child_id] = child_id +# Interpolates values from volume Gauss nodes to face nodes on one element. +#! format: off +@inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector, + A::TensorProductGaussFaceOperator{3, Projection{ApplyFaceWeights}}, + x::AbstractVector) where {ApplyFaceWeights} +#! format: on + @unpack interp_matrix_gauss_to_face_1d, face_indices_tensor_product = A + @unpack inv_volume_weights_1d, nnodes_1d, nfaces = A + + fill!(out_vec, zero(eltype(out_vec))) + + # As of Julia 1.9, Base.ReshapedArray does not produce allocations when setting values. + # Thus, Base.ReshapedArray should be used if you are setting values in the array. + # `reshape` is fine if you are only accessing values. + # Note that, for 3D GaussSBP nodes, the indexing is first in y, then x, then z. + out = Base.ReshapedArray(out_vec, (nnodes_1d, nnodes_1d, nnodes_1d), ()) + + if ApplyFaceWeights == true + @turbo for i in eachindex(x) + x[i] = x[i] * A.face_weights[i] end end - return nothing -end + # interpolation in the y-direction + @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face + index_left = face_indices_tensor_product[1, i, j, 2] + index_right = face_indices_tensor_product[2, i, j, 2] + for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes + out[jj, i, j] = out[jj, i, j] + + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] + out[jj, i, j] = out[jj, i, j] + + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] + end + end -# Refine given cells without rebalancing tree. -# -# Note: After a call to this method the tree may be unbalanced! -function refine_unbalanced!(t::AbstractTree, cell_ids, - sorted_unique_cell_ids = sort(unique(cell_ids))) - # Store actual ids refined cells (shifted due to previous insertions) - refined = zeros(Int, length(cell_ids)) - - # Loop over all cells that are to be refined - for (count, original_cell_id) in enumerate(sorted_unique_cell_ids) - # Determine actual cell id, taking into account previously inserted cells - n_children = n_children_per_cell(t) - cell_id = original_cell_id + (count - 1) * n_children - refined[count] = cell_id - - @assert !has_children(t, cell_id) "Non-leaf cell $cell_id cannot be refined" - - # Insert new cells directly behind parent (depth-first) - insert!(t, cell_id + 1, n_children) - - # Flip sign of refined cell such that we can easily find it later - t.original_cell_ids[cell_id] = -t.original_cell_ids[cell_id] - - # Initialize child cells (except neighbors) - for child in 1:n_children - child_id = cell_id + child - init_child!(t, cell_id, child, child_id) + # interpolation in the x-direction + @turbo for j in Base.OneTo(nnodes_1d), i in Base.OneTo(nnodes_1d) # loop over nodes in a face + index_left = face_indices_tensor_product[1, i, j, 1] + index_right = face_indices_tensor_product[2, i, j, 1] + for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes + out[i, jj, j] = out[i, jj, j] + + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] + out[i, jj, j] = out[i, jj, j] + + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] end + end - # Initialize child cells (only neighbors) - # This separate loop is required since init_child_neighbors requires initialized parent-child - # relationships - for child in 1:n_children - child_id = cell_id + child - init_child_neighbors!(t, cell_id, child, child_id) + # interpolation in the z-direction + @turbo for i in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d) # loop over nodes in a face + index_left = face_indices_tensor_product[1, i, j, 3] + index_right = face_indices_tensor_product[2, i, j, 3] + for jj in Base.OneTo(nnodes_1d) # loop over "line" of volume nodes + # The ordering (i,j) -> (j,i) needs to be reversed for this last face. + # This is due to way we define face nodes for Hex() types in StartUpDG.jl. + out[j, i, jj] = out[j, i, jj] + + interp_matrix_gauss_to_face_1d[1, jj] * x[index_left] + out[j, i, jj] = out[j, i, jj] + + interp_matrix_gauss_to_face_1d[2, jj] * x[index_right] end end - return refined + # apply inv(M) + @turbo for k in Base.OneTo(nnodes_1d), j in Base.OneTo(nnodes_1d), + i in Base.OneTo(nnodes_1d) + + out[i, j, k] = out[i, j, k] * inv_volume_weights_1d[i] * + inv_volume_weights_1d[j] * inv_volume_weights_1d[k] + end end -# Refine entire tree by one level -function refine!(t::AbstractTree) - cells = @trixi_timeit timer() "collect all leaf cells" leaf_cells(t) - @trixi_timeit timer() "refine!" refine!(t, cells, cells) +# For now, this is mostly the same as `create_cache` for DGMultiFluxDiff{<:Polynomial}. +# In the future, we may modify it so that we can specialize additional parts of GaussSBP() solvers. +function create_cache(mesh::DGMultiMesh, equations, + dg::DGMultiFluxDiff{<:GaussSBP, <:Union{Quad, Hex}}, RealT, + uEltype) + + # call general Polynomial flux differencing constructor + cache = invoke(create_cache, + Tuple{typeof(mesh), typeof(equations), + DGMultiFluxDiff, typeof(RealT), typeof(uEltype)}, + mesh, equations, dg, RealT, uEltype) + + rd = dg.basis + @unpack md = mesh + + # for change of basis prior to the volume integral and entropy projection + r1D, _ = StartUpDG.gauss_lobatto_quad(0, 0, polydeg(dg)) + rq1D, _ = StartUpDG.gauss_quad(0, 0, polydeg(dg)) + interp_matrix_lobatto_to_gauss_1D = polynomial_interpolation_matrix(r1D, rq1D) + interp_matrix_gauss_to_lobatto_1D = polynomial_interpolation_matrix(rq1D, r1D) + NDIMS = ndims(rd.element_type) + interp_matrix_lobatto_to_gauss = SimpleKronecker(NDIMS, + interp_matrix_lobatto_to_gauss_1D, + uEltype) + interp_matrix_gauss_to_lobatto = SimpleKronecker(NDIMS, + interp_matrix_gauss_to_lobatto_1D, + uEltype) + inv_gauss_weights = inv.(rd.wq) + + # specialized operators to perform tensor product interpolation to faces for Gauss nodes + interp_matrix_gauss_to_face = TensorProductGaussFaceOperator(Interpolation(), dg) + projection_matrix_gauss_to_face = TensorProductGaussFaceOperator(Projection{Static.False()}(), + dg) + + # `LIFT` matrix for Gauss nodes - this is equivalent to `projection_matrix_gauss_to_face` scaled by `diagm(rd.wf)`, + # where `rd.wf` are Gauss node face quadrature weights. + gauss_LIFT = TensorProductGaussFaceOperator(Projection{Static.True()}(), dg) + + nvars = nvariables(equations) + rhs_volume_local_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) + for _ in 1:Threads.nthreads()] + gauss_volume_local_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) + for _ in 1:Threads.nthreads()] + + return (; cache..., projection_matrix_gauss_to_face, gauss_LIFT, inv_gauss_weights, + rhs_volume_local_threaded, gauss_volume_local_threaded, + interp_matrix_lobatto_to_gauss, interp_matrix_gauss_to_lobatto, + interp_matrix_gauss_to_face, + create_cache(mesh, equations, dg.volume_integral, dg, RealT, uEltype)...) # add cache specialized on the volume integral end -# Refine given cells and rebalance tree. -# -# Note 1: Rebalancing is iterative, i.e., neighboring cells are refined if -# otherwise the 2:1 rule would be violated, which can cause more -# refinements. -# Note 2: Rebalancing currently only considers *Cartesian* neighbors, not diagonal neighbors! -function refine!(t::AbstractTree, cell_ids, - sorted_unique_cell_ids = sort(unique(cell_ids))) - # Reset original cell ids such that each cell knows its current id - reset_original_cell_ids!(t) - - # Refine all requested cells - refined = @trixi_timeit timer() "refine_unbalanced!" refine_unbalanced!(t, cell_ids, - sorted_unique_cell_ids) - refinement_count = length(refined) - - # Iteratively rebalance the tree until it does not change anymore - while length(refined) > 0 - refined = @trixi_timeit timer() "rebalance!" rebalance!(t, refined) - refinement_count += length(refined) +# by default, return an empty tuple for volume integral caches +create_cache(mesh, equations, volume_integral, dg, RealT, uEltype) = NamedTuple() + +# TODO: DGMulti. Address hard-coding of `entropy2cons!` and `cons2entropy!` for this function. +function entropy_projection!(cache, u, mesh::DGMultiMesh, equations, + dg::DGMultiFluxDiff{<:GaussSBP}) + rd = dg.basis + @unpack Vq = rd + @unpack VhP, entropy_var_values, u_values = cache + @unpack projected_entropy_var_values, entropy_projected_u_values = cache + @unpack interp_matrix_lobatto_to_gauss, interp_matrix_gauss_to_face = cache + + @threaded for e in eachelement(mesh, dg, cache) + apply_to_each_field(mul_by!(interp_matrix_lobatto_to_gauss), + view(u_values, :, e), view(u, :, e)) end - # Determine list of *original* cell ids that were refined - # Note: original_cell_ids contains the cell_id *before* refinement. At - # refinement, the refined cell's original_cell_ids value has its sign flipped - # to easily find it now. - refined_original_cells = @views(-t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0]) + # transform quadrature values to entropy variables + cons2entropy!(entropy_var_values, u_values, equations) - # Check if count of refinement cells matches information in original_cell_ids - @assert refinement_count==length(refined_original_cells) ("Mismatch in number of refined cells") + volume_indices = Base.OneTo(rd.Nq) + face_indices = (rd.Nq + 1):(rd.Nq + rd.Nfq) - return refined_original_cells -end - -# Refine all leaf cells with coordinates in a given rectangular box -function refine_box!(t::AbstractTree{NDIMS}, coordinates_min, - coordinates_max) where {NDIMS} - for dim in 1:NDIMS - @assert coordinates_min[dim] cell_coordinates(t, cell_id))) + # directly copy over volume values (no entropy projection required) + entropy_projected_volume_values = view(entropy_projected_u_values, volume_indices, + :) + @threaded for i in eachindex(u_values) + entropy_projected_volume_values[i] = u_values[i] end - # Refine cells - refine!(t, cells) -end + # transform entropy to conservative variables on face values + entropy_projected_face_values = view(entropy_projected_u_values, face_indices, :) + entropy2cons!(entropy_projected_face_values, entropy_var_face_values, equations) -# Convenience method for 1D -function refine_box!(t::AbstractTree{1}, coordinates_min::Real, coordinates_max::Real) - return refine_box!(t, [convert(Float64, coordinates_min)], - [convert(Float64, coordinates_max)]) + return nothing end -# Refine all leaf cells with coordinates in a given sphere -function refine_sphere!(t::AbstractTree{NDIMS}, center::SVector{NDIMS}, - radius) where {NDIMS} - @assert radius>=0 "Radius must be positive." +# Assumes cache.flux_face_values is already computed. +# Enables tensor product evaluation of `LIFT isa TensorProductGaussFaceOperator`. +function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations, + surface_integral::SurfaceIntegralWeakForm, + dg::DGMultiFluxDiff{<:GaussSBP}, cache) + (; gauss_LIFT, gauss_volume_local_threaded) = cache - # Find all leaf cells within sphere - cells = filter_leaf_cells(t) do cell_id - return sum(abs2, cell_coordinates(t, cell_id) - center) < radius^2 - end + @threaded for e in eachelement(mesh, dg, cache) - # Refine cells - refine!(t, cells) -end + # applies LIFT matrix, output is stored at Gauss nodes + gauss_volume_local = gauss_volume_local_threaded[Threads.threadid()] + apply_to_each_field(mul_by!(gauss_LIFT), gauss_volume_local, + view(cache.flux_face_values, :, e)) -# Convenience function to allow passing center as a tuple -function refine_sphere!(t::AbstractTree{NDIMS}, center::NTuple{NDIMS}, - radius) where {NDIMS} - refine_sphere!(t, SVector(center), radius) -end - -# For the given cell ids, check if neighbors need to be refined to restore a rebalanced tree. -# -# Note 1: Rebalancing currently only considers *Cartesian* neighbors, not diagonal neighbors! -# Note 2: The current algorithm assumes that a previous refinement step has -# created level differences of at most 2. That is, before the previous -# refinement step, the tree was balanced. -function rebalance!(t::AbstractTree, refined_cell_ids) - # Create buffer for newly refined cells - to_refine = zeros(Int, n_directions(t) * length(refined_cell_ids)) - count = 0 - - # Iterate over cell ids that have previously been refined - for cell_id in refined_cell_ids - # Go over all potential neighbors of child cell - for direction in eachdirection(t) - # Continue if refined cell has a neighbor in that direction - if has_neighbor(t, cell_id, direction) - continue - end - - # Continue if refined cell has no coarse neighbor, since that would - # mean it there is no neighbor in that direction at all (domain - # boundary) - if !has_coarse_neighbor(t, cell_id, direction) - continue - end - - # Otherwise, the coarse neighbor exists and is not refined, thus it must - # be marked for refinement - coarse_neighbor_id = t.neighbor_ids[direction, t.parent_ids[cell_id]] - count += 1 - to_refine[count] = coarse_neighbor_id + for i in eachindex(gauss_volume_local) + du[i, e] = du[i, e] + gauss_volume_local[i] end end - - # Finally, refine all marked cells... - refined = refine_unbalanced!(t, unique(to_refine[1:count])) - - # ...and return list of refined cells - return refined end -# Refine given cells without rebalancing tree. -# -# Note: After a call to this method the tree may be unbalanced! -# function refine_unbalanced!(t::AbstractTree, cell_ids) end - -# Wrap single-cell refinements such that `sort(...)` does not complain -refine_unbalanced!(t::AbstractTree, cell_id::Int) = refine_unbalanced!(t, [cell_id]) - -# Coarsen entire tree by one level -function coarsen!(t::AbstractTree) - # Special case: if there is only one cell (root), there is nothing to do - if length(t) == 1 - return Int[] +@inline function flux_differencing_kernel!(du, u, element, mesh::DGMultiMesh, + have_nonconservative_terms, equations, + volume_flux, dg::DGMultiFluxDiff{<:GaussSBP}, + cache, alpha = true) + fluxdiff_local = cache.fluxdiff_local_threaded[Threads.threadid()] + fill!(fluxdiff_local, zero(eltype(fluxdiff_local))) + u_local = view(cache.entropy_projected_u_values, :, element) + + local_flux_differencing!(fluxdiff_local, u_local, element, + have_nonconservative_terms, + volume_flux, has_sparse_operators(dg), + mesh, equations, dg, cache) + + # convert `fluxdiff_local::Vector{<:SVector}` to `rhs_local::StructArray{<:SVector}` + # for faster performance when using `apply_to_each_field`. + rhs_local = cache.rhs_local_threaded[Threads.threadid()] + for i in Base.OneTo(length(fluxdiff_local)) + rhs_local[i] = fluxdiff_local[i] end - # Get list of unique parent ids for all leaf cells - parent_ids = unique(t.parent_ids[leaf_cells(t)]) - coarsen!(t, parent_ids) + project_rhs_to_gauss_nodes!(du, rhs_local, element, mesh, dg, cache, alpha) end -# Coarsen given *parent* cells (= these cells must have children who are all -# leaf cells) while retaining a balanced tree. -# -# A cell to be coarsened might cause an unbalanced tree if the neighboring cell -# was already refined. Since it is generally not desired that cells are -# coarsened without specifically asking for it, these cells will then *not* be -# coarsened. -function coarsen!(t::AbstractTree, cell_ids::AbstractArray{Int}) - # Return early if array is empty - if length(cell_ids) == 0 - return Int[] +function project_rhs_to_gauss_nodes!(du, rhs_local, element, mesh::DGMultiMesh, + dg::DGMulti, cache, alpha = true) + + # Here, we exploit that under a Gauss nodal basis the structure of the projection + # matrix `Ph = [diagm(1 ./ wq), projection_matrix_gauss_to_face]` such that + # `Ph * [u; uf] = (u ./ wq) + projection_matrix_gauss_to_face * uf`. + volume_indices = Base.OneTo(dg.basis.Nq) + face_indices = (dg.basis.Nq + 1):(dg.basis.Nq + dg.basis.Nfq) + local_volume_flux = view(rhs_local, volume_indices) + local_face_flux = view(rhs_local, face_indices) + + # initialize rhs_volume_local = projection_matrix_gauss_to_face * local_face_flux + rhs_volume_local = cache.rhs_volume_local_threaded[Threads.threadid()] + apply_to_each_field(mul_by!(cache.projection_matrix_gauss_to_face), + rhs_volume_local, local_face_flux) + + # accumulate volume contributions at Gauss nodes + for i in eachindex(rhs_volume_local) + du_local = rhs_volume_local[i] + + local_volume_flux[i] * cache.inv_gauss_weights[i] + du[i, element] = du[i, element] + alpha * du_local end +end - # Reset original cell ids such that each cell knows its current id - reset_original_cell_ids!(t) - - # To maximize the number of cells that may be coarsened, start with the cells at the highest level - sorted_by_level = sort(cell_ids, by = i -> t.levels[i]) - - # Keep track of number of cells that were actually coarsened - n_coarsened = 0 - - # Local function to adjust cell ids after some cells have been removed - function adjust_cell_ids!(cell_ids, coarsened_cell_id, count) - for (id, cell_id) in enumerate(cell_ids) - if cell_id > coarsened_cell_id - cell_ids[id] = cell_id - count - end - end +function calc_volume_integral!(du, u, mesh::DGMultiMesh, + have_nonconservative_terms, equations, + volume_integral::VolumeIntegralFluxDifferencing, + dg::DGMultiFluxDiff{<:GaussSBP}, cache) + @threaded for e in eachelement(mesh, dg, cache) + flux_differencing_kernel!(du, u, e, mesh, + have_nonconservative_terms, equations, + volume_integral.volume_flux, dg, cache) end +end - # Iterate backwards over cells to coarsen - while true - # Retrieve next cell or quit - if length(sorted_by_level) > 0 - coarse_cell_id = pop!(sorted_by_level) - else - break - end - - # Ensure that cell has children (violation is an error) - if !has_children(t, coarse_cell_id) - error("cell is leaf and cannot be coarsened to: $coarse_cell_id") - end +# interpolate back to Lobatto nodes after applying the inverse Jacobian at Gauss points +function invert_jacobian_and_interpolate!(du, mesh::DGMultiMesh, equations, + dg::DGMultiFluxDiff{<:GaussSBP}, cache; + scaling = -1) + (; interp_matrix_gauss_to_lobatto, rhs_volume_local_threaded, invJ) = cache - # Ensure that all child cells are leaf cells (violation is an error) - for child in 1:n_children_per_cell(t) - if has_child(t, coarse_cell_id, child) - if !is_leaf(t, t.child_ids[child, coarse_cell_id]) - error("cell $coarse_cell_id has child cell at position $child that is not a leaf cell") - end - end - end + @threaded for e in eachelement(mesh, dg, cache) + rhs_volume_local = rhs_volume_local_threaded[Threads.threadid()] - # Check if coarse cell has refined neighbors that would prevent coarsening - skip = false - # Iterate over all children (which are to be removed) - for child in 1:n_children_per_cell(t) - # Continue if child does not exist - if !has_child(t, coarse_cell_id, child) - continue - end - child_id = t.child_ids[child, coarse_cell_id] - - # Go over all neighbors of child cell. If it has a neighbor that is *not* - # a sibling and that is not a leaf cell, we cannot coarsen its parent - # without creating an unbalanced tree. - for direction in eachdirection(t) - # Continue if neighbor would be a sibling - if has_sibling(child, direction) - continue - end - - # Continue if child cell has no neighbor in that direction - if !has_neighbor(t, child_id, direction) - continue - end - neighbor_id = t.neighbor_ids[direction, child_id] - - if !has_children(t, neighbor_id) - continue - end - - # If neighbor is not a sibling, is existing, and has children, do not coarsen - skip = true - break - end - end - # Skip if a neighboring cell prevents coarsening - if skip - continue + # At this point, `rhs_volume_local` should still be stored at Gauss points. + # We scale it by the inverse Jacobian before transforming back to Lobatto. + for i in eachindex(rhs_volume_local) + rhs_volume_local[i] = du[i, e] * invJ[i, e] * scaling end - # Flip sign of cell to be coarsened to such that we can easily find it - t.original_cell_ids[coarse_cell_id] = -t.original_cell_ids[coarse_cell_id] - - # If a coarse cell has children that are all leaf cells, they must follow - # immediately due to depth-first ordering of the tree - count = n_children(t, coarse_cell_id) - @assert count==n_children_per_cell(t) "cell $coarse_cell_id does not have all child cells" - remove_shift!(t, coarse_cell_id + 1, coarse_cell_id + count) - - # Take into account shifts in tree that alters cell ids - adjust_cell_ids!(sorted_by_level, coarse_cell_id, count) - - # Keep track of number of coarsened cells - n_coarsened += 1 + # Interpolate result back to Lobatto nodes for ease of analysis, visualization + apply_to_each_field(mul_by!(interp_matrix_gauss_to_lobatto), + view(du, :, e), rhs_volume_local) end - - # Determine list of *original* cell ids that were coarsened to - # Note: original_cell_ids contains the cell_id *before* coarsening. At - # coarsening, the coarsened parent cell's original_cell_ids value has its sign flipped - # to easily find it now. - @views coarsened_original_cells = (-t.original_cell_ids[1:length(t)][t.original_cell_ids[1:length(t)] .< 0]) - - # Check if count of coarsened cells matches information in original_cell_ids - @assert n_coarsened==length(coarsened_original_cells) ("Mismatch in number of coarsened cells") - - return coarsened_original_cells end -# Wrap single-cell coarsening such that `sort(...)` does not complain -coarsen!(t::AbstractTree, cell_id::Int) = coarsen!(t::AbstractTree, [cell_id]) - -# Coarsen all viable parent cells with coordinates in a given rectangular box -function coarsen_box!(t::AbstractTree{NDIMS}, coordinates_min::AbstractArray{Float64}, - coordinates_max::AbstractArray{Float64}) where {NDIMS} - for dim in 1:NDIMS - @assert coordinates_min[dim] cell_coordinates(t, cell_id))) + # `du` is stored at Gauss nodes here + @trixi_timeit timer() "volume integral" begin + calc_volume_integral!(du, u, mesh, + have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache) end - # Get list of unique parent ids for all leaf cells - parent_ids = unique(t.parent_ids[leaves]) - - # Filter parent ids to be within box - parents = filter(parent_ids) do cell_id - return (all(coordinates_min .< cell_coordinates(t, cell_id)) && - all(coordinates_max .> cell_coordinates(t, cell_id))) + # the following functions are the same as in VolumeIntegralWeakForm, and can be reused from dg.jl + @trixi_timeit timer() "interface flux" begin + calc_interface_flux!(cache, dg.surface_integral, mesh, + have_nonconservative_terms(equations), equations, dg) end - # Coarsen cells - coarsen!(t, parents) -end - -# Convenience method for 1D -function coarsen_box!(t::AbstractTree{1}, coordinates_min::Real, coordinates_max::Real) - return coarsen_box!(t, [convert(Float64, coordinates_min)], - [convert(Float64, coordinates_max)]) -end - -# Return coordinates of a child cell based on its relative position to the parent. -function child_coordinates(::AbstractTree{NDIMS}, parent_coordinates, - parent_length::Number, child::Int) where {NDIMS} - # Calculate length of child cells - child_length = parent_length / 2 - return SVector(ntuple(d -> parent_coordinates[d] + - child_sign(child, d) * child_length / 2, Val(NDIMS))) -end - -# Reset range of cells to values that are prone to cause errors as soon as they are used. -# -# Rationale: If an invalid cell is accidentally used, we want to know it as soon as possible. -# function invalidate!(t::AbstractTree, first::Int, last::Int) end -invalidate!(t::AbstractTree, id::Int) = invalidate!(t, id, id) -invalidate!(t::AbstractTree) = invalidate!(t, 1, length(t)) - -# Delete connectivity with parents/children/neighbors before cells are erased -function delete_connectivity!(t::AbstractTree, first::Int, last::Int) - @assert first > 0 - @assert first <= last - @assert last <= t.capacity + 1 - - # Iterate over all cells - for cell_id in first:last - # Delete connectivity from parent cell - if has_parent(t, cell_id) - parent_id = t.parent_ids[cell_id] - for child in 1:n_children_per_cell(t) - if t.child_ids[child, parent_id] == cell_id - t.child_ids[child, parent_id] = 0 - break - end - end - end - - # Delete connectivity from child cells - for child in 1:n_children_per_cell(t) - if has_child(t, cell_id, child) - t.parent_ids[t._child_ids[child, cell_id]] = 0 - end - end - - # Delete connectivity from neighboring cells - for direction in eachdirection(t) - if has_neighbor(t, cell_id, direction) - t.neighbor_ids[opposite_direction(direction), t.neighbor_ids[direction, cell_id]] = 0 - end - end + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux!(cache, t, boundary_conditions, mesh, + have_nonconservative_terms(equations), equations, dg) end -end - -# Move connectivity with parents/children/neighbors after cells have been moved -function move_connectivity!(t::AbstractTree, first::Int, last::Int, destination::Int) - @assert first > 0 - @assert first <= last - @assert last <= t.capacity + 1 - @assert destination > 0 - @assert destination <= t.capacity + 1 - - # Strategy - # 1) Loop over moved cells (at target location) - # 2) Check if parent/children/neighbors connections are to a cell that was moved - # a) if cell was moved: apply offset to current cell - # b) if cell was not moved: go to connected cell and update connectivity there - - offset = destination - first - has_moved(n) = (first <= n <= last) - - for source in first:last - target = source + offset - - # Update parent - if has_parent(t, target) - # Get parent cell - parent_id = t.parent_ids[target] - if has_moved(parent_id) - # If parent itself was moved, just update parent id accordingly - t.parent_ids[target] += offset - else - # If parent was not moved, update its corresponding child id - for child in 1:n_children_per_cell(t) - if t.child_ids[child, parent_id] == source - t.child_ids[child, parent_id] = target - end - end - end - end - - # Update children - for child in 1:n_children_per_cell(t) - if has_child(t, target, child) - # Get child cell - child_id = t.child_ids[child, target] - if has_moved(child_id) - # If child itself was moved, just update child id accordingly - t.child_ids[child, target] += offset - else - # If child was not moved, update its parent id - t.parent_ids[child_id] = target - end - end - end - # Update neighbors - for direction in eachdirection(t) - if has_neighbor(t, target, direction) - # Get neighbor cell - neighbor_id = t.neighbor_ids[direction, target] - if has_moved(neighbor_id) - # If neighbor itself was moved, just update neighbor id accordingly - t.neighbor_ids[direction, target] += offset - else - # If neighbor was not moved, update its opposing neighbor id - t.neighbor_ids[opposite_direction(direction), neighbor_id] = target - end - end - end + # `du` is stored at Gauss nodes here + @trixi_timeit timer() "surface integral" begin + calc_surface_integral!(du, u, mesh, equations, + dg.surface_integral, dg, cache) end -end -# Raw copy operation for ranges of cells. -# -# This method is used by the higher-level copy operations for AbstractContainer -# function raw_copy!(target::AbstractTree, source::AbstractTree, first::Int, last::Int, destination::Int) end + # invert Jacobian and map `du` from Gauss to Lobatto nodes + @trixi_timeit timer() "Jacobian" begin + invert_jacobian_and_interpolate!(du, mesh, equations, dg, cache) + end -# Reset data structures by recreating all internal storage containers and invalidating all elements -# function reset_data_structures!(t::AbstractTree{NDIMS}) where NDIMS end + @trixi_timeit timer() "source terms" begin + calc_sources!(du, u, t, source_terms, mesh, equations, dg, cache) + end + return nothing +end end # @muladd From 3cc6140a026771788aee5e58692bc94410c58e2a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 8 Apr 2024 09:32:27 +0200 Subject: [PATCH 5/9] revert --- src/solvers/dgsem_tree/dg_2d_parallel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index 746a9802d66..157d462aa2f 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -60,7 +60,7 @@ end function start_mpi_send!(mpi_cache::MPICache, mesh, equations, dg, cache) data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1) - for d in eachindex(mpi_cache.mpi_neighbor_ranks) + for d in 1:length(mpi_cache.mpi_neighbor_ranks) send_buffer = mpi_cache.mpi_send_buffers[d] for (index, interface) in enumerate(mpi_cache.mpi_neighbor_interfaces[d]) From d2350854f3924841ee0f4f5e9ec9a55c8ed5d0d0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 8 Apr 2024 09:36:21 +0200 Subject: [PATCH 6/9] revert --- ext/TrixiMakieExt.jl | 2 +- src/semidiscretization/semidiscretization_coupled.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ext/TrixiMakieExt.jl b/ext/TrixiMakieExt.jl index 1684ea2dc97..301a7656da9 100644 --- a/ext/TrixiMakieExt.jl +++ b/ext/TrixiMakieExt.jl @@ -187,7 +187,7 @@ function iplot(pd::PlotData2DTriangulated; fig = Makie.Figure() # Set up options for the drop-down menu - menu_options = [zip(variable_names, eachindex(variable_names))...] + menu_options = [zip(variable_names, 1:length(variable_names))...] menu = Makie.Menu(fig, options = menu_options) # Initialize toggle switches for viewing the mesh diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 9ed01644a52..dc21dbe9a1e 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -32,14 +32,14 @@ function SemidiscretizationCoupled(semis...) # Number of coefficients for each semidiscretization n_coefficients = zeros(Int, length(semis)) - for i in eachindex(semis) + for i in 1:length(semis) _, equations, _, _ = mesh_equations_solver_cache(semis[i]) n_coefficients[i] = ndofs(semis[i]) * nvariables(equations) end # Compute range of coefficients associated with each semidiscretization and allocate coupled BCs u_indices = Vector{UnitRange{Int}}(undef, length(semis)) - for i in eachindex(semis) + for i in 1:length(semis) offset = sum(n_coefficients[1:(i - 1)]) + 1 u_indices[i] = range(offset, length = n_coefficients[i]) From 749e55f5dedf9eae6ba3858e0e6a576c78c6c05e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 8 Apr 2024 09:41:52 +0200 Subject: [PATCH 7/9] fmt --- src/solvers/dgmulti/flux_differencing_gauss_sbp.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index 7c9227b5c3c..892db073e06 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -175,8 +175,8 @@ end # Interpolates values from volume Gauss nodes to face nodes on one element. #! format: off @inline function tensor_product_gauss_face_operator!(out::AbstractVector, - A::TensorProductGaussFaceOperator{3, Interpolation}, - x::AbstractVector) + A::TensorProductGaussFaceOperator{3, Interpolation}, + x::AbstractVector) #! format: on (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A (; nnodes_1d) = A @@ -228,8 +228,8 @@ end # Projects face node values to volume Gauss nodes on one element. #! format: off @inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector, - A::TensorProductGaussFaceOperator{2, Projection{ApplyFaceWeights}}, - x::AbstractVector) where {ApplyFaceWeights} + A::TensorProductGaussFaceOperator{2, Projection{ApplyFaceWeights}}, + x::AbstractVector) where {ApplyFaceWeights} #! format: on (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A (; inv_volume_weights_1d, nnodes_1d) = A @@ -282,8 +282,8 @@ end # Interpolates values from volume Gauss nodes to face nodes on one element. #! format: off @inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector, - A::TensorProductGaussFaceOperator{3, Projection{ApplyFaceWeights}}, - x::AbstractVector) where {ApplyFaceWeights} + A::TensorProductGaussFaceOperator{3, Projection{ApplyFaceWeights}}, + x::AbstractVector) where {ApplyFaceWeights} #! format: on @unpack interp_matrix_gauss_to_face_1d, face_indices_tensor_product = A @unpack inv_volume_weights_1d, nnodes_1d, nfaces = A From 770ac60176fe0f826585185b15d7da3a6693c23b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 8 Apr 2024 15:01:27 +0200 Subject: [PATCH 8/9] fmt --- src/solvers/dgmulti/flux_differencing_gauss_sbp.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index 892db073e06..9059caf87f6 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -136,8 +136,8 @@ end #! format: off # Interpolates values from volume Gauss nodes to face nodes on one element. @inline function tensor_product_gauss_face_operator!(out::AbstractVector, - A::TensorProductGaussFaceOperator{2, Interpolation}, - x_in::AbstractVector) + A::TensorProductGaussFaceOperator{2, Interpolation}, + x_in::AbstractVector) #! format: on (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A (; nnodes_1d) = A From acac1b748a263a3366bdce3ca6627093a6c1cca8 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 9 Apr 2024 09:05:56 +0200 Subject: [PATCH 9/9] revert --- src/callbacks_step/time_series_dg_unstructured.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/time_series_dg_unstructured.jl b/src/callbacks_step/time_series_dg_unstructured.jl index e9a18fcd3bb..85427f1273a 100644 --- a/src/callbacks_step/time_series_dg_unstructured.jl +++ b/src/callbacks_step/time_series_dg_unstructured.jl @@ -118,7 +118,7 @@ function calc_minimum_surface_distance(point, node_coordinates, n = nnodes(dg) min_distance2 = Inf * ones(eltype(mesh.corners), length(mesh)) indices = zeros(Int, length(mesh), 2) - for k in eachindex(mesh) + for k in 1:length(mesh) # used to ensure that only boundary points are used on_surface = MVector(false, false) for j in 1:n