Implement suggestions

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

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

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
 

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

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)

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)

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,