From e71913fb37fa7493c41c39c497a9287d51372ba4 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 9 Jun 2021 18:08:31 +0200 Subject: [PATCH] Implement suggestions --- src/auxiliary/auxiliary.jl | 2 ++ src/callbacks_step/amr.jl | 19 +++++++++----- src/callbacks_step/analysis.jl | 2 +- src/mesh/p4est_mesh.jl | 22 +++++++++++------ src/solvers/dg_p4est/containers.jl | 15 ++++++------ src/solvers/dg_p4est/containers_2d.jl | 30 +++++++++++------------ src/solvers/dg_unstructured_quad/dg_2d.jl | 6 ----- 7 files changed, 53 insertions(+), 43 deletions(-) diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 8318ba9d6bd..65f8a2f3d33 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -228,4 +228,6 @@ function init_p4est() # Initialize p4est with log level ERROR to prevent a lot of output in AMR simulations p4est_init(C_NULL, SC_LP_ERROR) + + return nothing end diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 55ea59bfe93..156c2b32dd7 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -169,7 +169,7 @@ end kwargs...) # Note that we don't `wrap_array` the vector `u_ode` to be able to `resize!` # it when doing AMR while still dispatching on the `mesh` etc. - amr_callback(u_ode, mesh_equations_solver_cache(semi)..., t, iter; kwargs...) + amr_callback(u_ode, mesh_equations_solver_cache(semi)..., semi, t, iter; kwargs...) end @@ -181,7 +181,7 @@ end # `passive_args` is expected to be an iterable of `Tuple`s of the form # `(p_u_ode, p_mesh, p_equations, p_dg, p_cache)`. function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, - equations, dg::DG, cache, + equations, dg::DG, cache, semi, t, iter; only_refine=false, only_coarsen=false, passive_args=()) @@ -322,14 +322,13 @@ end # Copy controller values to quad user data storage, will be called below function copy_to_quad_iter_volume(info, user_data) # Load tree from global trees array, one-based indexing - tree = load_sc_array(p4est_tree_t, info.p4est.trees, info.treeid + 1) + tree = unsafe_load_sc(p4est_tree_t, info.p4est.trees, info.treeid + 1) # Quadrant numbering offset of this quadrant offset = tree.quadrants_offset # Global quad ID quad_id = offset + info.quadid # Access user_data = lambda - # TODO P4EST Is lambda always a Vector{Int}? user_data_ptr = Ptr{Int}(user_data) # Load controller_value = lambda[quad_id + 1] controller_value = unsafe_load(user_data_ptr, quad_id + 1) @@ -343,7 +342,7 @@ function copy_to_quad_iter_volume(info, user_data) end function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh, - equations, dg::DG, cache, + equations, dg::DG, cache, semi, t, iter; only_refine=false, only_coarsen=false, passive_args=()) @@ -362,6 +361,9 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh, # Copy controller value of each quad to the quad's user data storage iter_volume_c = @cfunction(copy_to_quad_iter_volume, Cvoid, (Ptr{p4est_iter_volume_info_t}, Ptr{Cvoid})) + # The pointer to lambda will be interpreted as Ptr{Int} above + @assert lambda isa Vector{Int} + GC.@preserve lambda begin p4est_iterate(mesh.p4est, C_NULL, # ghost layer @@ -408,6 +410,11 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh, if has_changed # TODO: Taal decide, where shall we set this? # don't set it to has_changed since there can be changes from earlier calls mesh.unsaved_changes = true + + if semi.boundary_conditions isa UnstructuredQuadSortedBoundaryTypes + # Reinitialize boundary types container because boundaries may have changed. + initialize!(boundary_conditions, cache) + end end # Return true if there were any cells coarsened or refined, otherwise false @@ -531,7 +538,7 @@ end function extract_levels_iter_volume(info, user_data) # Load tree from global trees array, one-based indexing - tree = load_sc_array(p4est_tree_t, info.p4est.trees, info.treeid + 1) + tree = unsafe_load_sc(p4est_tree_t, info.p4est.trees, info.treeid + 1) # Quadrant numbering offset of this quadrant offset = tree.quadrants_offset # Global quad ID diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 80c95d28e73..c306d965c67 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -385,7 +385,7 @@ function print_amr_information(callbacks, mesh::P4estMesh, solver, cache) elements_per_level = zeros(P4EST_MAXLEVEL + 1) - for tree in convert_sc_array(p4est_tree_t, mesh.p4est.trees) + for tree in unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees) elements_per_level .+= tree.quadrants_per_level end diff --git a/src/mesh/p4est_mesh.jl b/src/mesh/p4est_mesh.jl index e19bfa73513..42156cdbbc7 100644 --- a/src/mesh/p4est_mesh.jl +++ b/src/mesh/p4est_mesh.jl @@ -457,8 +457,12 @@ function coarsen_fn(p4est, which_tree, quadrants_ptr) quadrants = unsafe_wrap(Array, quadrants_ptr, 4) # Controller value has been copied to the quadrant's user data storage before. - # Unpack quadrant's user data ([global quad ID, controller_value]) - if all(i -> unsafe_load(Ptr{Int}(quadrants[i].p.user_data), 2) < 0, eachindex(quadrants)) + # Load controller value from quadrant's user data ([global quad ID, controller_value]). + controller_value(i) = unsafe_load(Ptr{Int}(quadrants[i].p.user_data), 2) + + # p4est calls this function for each 2^ndims quads that could be coarsened to a single one. + # Only coarsen if all these 2^ndims quads have been marked for coarsening. + if all(i -> controller_value(i) < 0, eachindex(quadrants)) # return true (coarsen) return Cint(1) else @@ -482,9 +486,11 @@ function coarsen!(mesh::P4estMesh) @timed timer() "coarsen!" p4est_coarsen(mesh.p4est, false, coarsen_fn_c, init_fn_c) + # IDs of newly created cells (one-based) new_cells = collect_new_cells(mesh) + # Old IDs of cells that have been coarsened (one-based) coarsened_cells_vec = collect_changed_cells(mesh, original_n_cells) - # 2^NDIMS changed cells should have been coarsened to one new cell. + # 2^ndims changed cells should have been coarsened to one new cell. # This matrix will store the IDs of all cells that have been coarsened to cell new_cells[i] # in the i-th column. coarsened_cells = reshape(coarsened_cells_vec, 2^ndims(mesh), length(new_cells)) @@ -520,8 +526,8 @@ end # Copy global quad ID to quad's user data storage, will be called below function save_original_id_iter_volume(info, user_data) - # Global trees array, one-based indexing - tree = load_sc_array(p4est_tree_t, info.p4est.trees, info.treeid + 1) + # Load tree from global trees array, one-based indexing + tree = unsafe_load_sc(p4est_tree_t, info.p4est.trees, info.treeid + 1) # Quadrant numbering offset of this quadrant offset = tree.quadrants_offset # Global quad ID @@ -551,7 +557,7 @@ end # Extract information about which cells have been changed function collect_changed_iter_volume(info, user_data) # The original element ID has been saved to user_data before. - # Unpack quadrant's user data ([global quad ID, controller_value]). + # Load original quad ID from quad's user data ([global quad ID, controller_value]). quad_data_ptr = Ptr{Int}(info.quad.p.user_data) original_id = unsafe_load(quad_data_ptr, 1) @@ -600,8 +606,8 @@ function collect_new_iter_volume(info, user_data) # original_id of cells that have been newly created is -1 if original_id < 0 - # Global trees array, one-based indexing - tree = load_sc_array(p4est_tree_t, info.p4est.trees, info.treeid + 1) + # Load tree from global trees array, one-based indexing + tree = unsafe_load_sc(p4est_tree_t, info.p4est.trees, info.treeid + 1) # Quadrant numbering offset of this quadrant offset = tree.quadrants_offset # Global quad ID diff --git a/src/solvers/dg_p4est/containers.jl b/src/solvers/dg_p4est/containers.jl index 341c4ebfafa..de24bb03569 100644 --- a/src/solvers/dg_p4est/containers.jl +++ b/src/solvers/dg_p4est/containers.jl @@ -339,8 +339,8 @@ end function count_interfaces_iter_face(info, user_data) if info.sides.elem_count == 2 # Extract interface data - sides = (load_sc_array(p4est_iter_face_side_t, info.sides, 1), - load_sc_array(p4est_iter_face_side_t, info.sides, 2)) + sides = (unsafe_load_sc(p4est_iter_face_side_t, info.sides, 1), + unsafe_load_sc(p4est_iter_face_side_t, info.sides, 2)) if sides[1].is_hanging == 0 && sides[2].is_hanging == 0 # No hanging nodes => normal interface @@ -370,8 +370,8 @@ end function count_mortars_iter_face(info, user_data) if info.sides.elem_count == 2 # Extract interface data - sides = (load_sc_array(p4est_iter_face_side_t, info.sides, 1), - load_sc_array(p4est_iter_face_side_t, info.sides, 2)) + sides = (unsafe_load_sc(p4est_iter_face_side_t, info.sides, 1), + unsafe_load_sc(p4est_iter_face_side_t, info.sides, 2)) if sides[1].is_hanging != 0 || sides[2].is_hanging != 0 # Hanging nodes => mortar @@ -431,8 +431,8 @@ function iterate_faces(mesh::P4estMesh, iter_face_c, user_data) end -# Convert sc_array to Julia array of the specified type -function convert_sc_array(::Type{T}, sc_array) where T +# Convert sc_array of type T to Julia array +function unsafe_wrap_sc(::Type{T}, sc_array) where T element_count = sc_array.elem_count element_size = sc_array.elem_size @@ -442,7 +442,8 @@ function convert_sc_array(::Type{T}, sc_array) where T end -function load_sc_array(::Type{T}, sc_array, i=1) where T +# Load the ith element (1-indexed) of an sc array of type T +function unsafe_load_sc(::Type{T}, sc_array, i=1) where T element_size = sc_array.elem_size @assert element_size == sizeof(T) diff --git a/src/solvers/dg_p4est/containers_2d.jl b/src/solvers/dg_p4est/containers_2d.jl index 232583bc884..23cc4111fda 100644 --- a/src/solvers/dg_p4est/containers_2d.jl +++ b/src/solvers/dg_p4est/containers_2d.jl @@ -31,11 +31,11 @@ function calc_node_coordinates!(node_coordinates, p4est_root_len = 1 << P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) - trees = convert_sc_array(p4est_tree_t, mesh.p4est.trees) + trees = unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees) for tree in eachindex(trees) offset = trees[tree].quadrants_offset - quadrants = convert_sc_array(p4est_quadrant_t, trees[tree].quadrants) + quadrants = unsafe_wrap_sc(p4est_quadrant_t, trees[tree].quadrants) for i in eachindex(quadrants) element = offset + i @@ -69,8 +69,8 @@ function init_interfaces_iter_face(info, user_data) return nothing end - sides = (load_sc_array(p4est_iter_face_side_t, info.sides, 1), - load_sc_array(p4est_iter_face_side_t, info.sides, 2)) + sides = (unsafe_load_sc(p4est_iter_face_side_t, info.sides, 1), + unsafe_load_sc(p4est_iter_face_side_t, info.sides, 2)) if sides[1].is_hanging == true || sides[2].is_hanging == true # Mortar, no normal interface @@ -91,9 +91,9 @@ end # Function barrier for type stability function init_interfaces_iter_face_inner(info, sides, interfaces, interface_id, mesh) - # Global trees array, one-based indexing - trees = (load_sc_array(p4est_tree_t, mesh.p4est.trees, sides[1].treeid + 1), - load_sc_array(p4est_tree_t, mesh.p4est.trees, sides[2].treeid + 1)) + # Load local trees from global trees array, one-based indexing + trees = (unsafe_load_sc(p4est_tree_t, mesh.p4est.trees, sides[1].treeid + 1), + unsafe_load_sc(p4est_tree_t, mesh.p4est.trees, sides[2].treeid + 1)) # Quadrant numbering offsets of the quadrants at this interface offsets = SVector(trees[1].quadrants_offset, trees[2].quadrants_offset) @@ -179,9 +179,9 @@ end # Function barrier for type stability function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh) # Extract boundary data - side = load_sc_array(p4est_iter_face_side_t, info.sides) - # Global trees array, one-based indexing - tree = load_sc_array(p4est_tree_t, mesh.p4est.trees, side.treeid + 1) + side = unsafe_load_sc(p4est_iter_face_side_t, info.sides) + # Load tree from global trees array, one-based indexing + tree = unsafe_load_sc(p4est_tree_t, mesh.p4est.trees, side.treeid + 1) # Quadrant numbering offset of this quadrant offset = tree.quadrants_offset @@ -238,8 +238,8 @@ function init_mortars_iter_face(info, user_data) return nothing end - sides = (load_sc_array(p4est_iter_face_side_t, info.sides, 1), - load_sc_array(p4est_iter_face_side_t, info.sides, 2)) + sides = (unsafe_load_sc(p4est_iter_face_side_t, info.sides, 1), + unsafe_load_sc(p4est_iter_face_side_t, info.sides, 2)) if sides[1].is_hanging == false && sides[2].is_hanging == false # Normal interface, no mortar @@ -260,9 +260,9 @@ end # Function barrier for type stability function init_mortars_iter_face_inner(info, sides, mortars, mortar_id, mesh) - # Global trees array, one-based indexing - trees = (load_sc_array(p4est_tree_t, mesh.p4est.trees, sides[1].treeid + 1), - load_sc_array(p4est_tree_t, mesh.p4est.trees, sides[2].treeid + 1)) + # Load local trees from global trees array, one-based indexing + trees = (unsafe_load_sc(p4est_tree_t, mesh.p4est.trees, sides[1].treeid + 1), + unsafe_load_sc(p4est_tree_t, mesh.p4est.trees, sides[2].treeid + 1)) # Quadrant numbering offsets of the quadrants at this interface offsets = SVector(trees[1].quadrants_offset, trees[2].quadrants_offset) diff --git a/src/solvers/dg_unstructured_quad/dg_2d.jl b/src/solvers/dg_unstructured_quad/dg_2d.jl index 77ce7afcd7b..697d2e6ac05 100644 --- a/src/solvers/dg_unstructured_quad/dg_2d.jl +++ b/src/solvers/dg_unstructured_quad/dg_2d.jl @@ -307,12 +307,6 @@ end function calc_boundary_flux!(cache, t, boundary_conditions, mesh::Union{UnstructuredQuadMesh, P4estMesh}, equations, dg::DG) - if nboundaries(cache.boundaries) != nboundaries(boundary_conditions) - # When using AMR, the boundary container may have been resized. - # In this case, the boundary types container must be reinitialized as well. - initialize!(boundary_conditions, cache) - end - @unpack boundary_condition_types, boundary_indices = boundary_conditions calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,