Error for 3D when using p4est 2.8.5

[JordiManyer]

Tests for 3D are failing very early on when using more modern p4est releases. Note: I am not sure for which version this starts failing, I've only tested versions 2.2 and 2.8.5.

So here is my analysis of the issue:

In p4est, the p8est_tree_t type is defined here, and it's variable quadrants_per_level is set as an array of length P8EST_MAXLEVEL+1.

In P4est_wrapper.jl, the structure that is meant to replicate it is defined here. The maximum number of levels is set to 19.

For p4est v2.8.5, we have that P8EST_MAXLEVEL=30, whereas in p4est 2.2 it was indeed 19. I suspect this is meant to increase consistency between the p4est and p8est implementations (now both have same number of max levels).

So a new version of P4est_wrapper is necessary to correct these changes.

[wei3li]

Hi @amartinhuertas and @JordiManyer,

I am not sure if the crash I encountered is relevant to this issue. When I was using the "hand-drawn" 3D coarse model to construct an OctreeDistributedDiscreteModel, the code crashed. Here is the code I was running. Note that the init_coarse_model is the same as setup_model in GridapP4est.jl/test/CoarseDiscreteModelsTools.jl

using Gridap, GridapP4est, GridapDistributed, PartitionedArrays
import MPI
import FillArrays: Fill

function init_coarse_model()
  ptr = [1, 9, 17]
  data = [1, 2, 3, 4, 5, 6, 7, 8, 2, 9, 4, 10, 6, 11, 8, 12]
  cell_vertex_lids = Gridap.Arrays.Table(data, ptr)
  node_coordinates = Vector{Point{3,Float64}}(undef, 12)
  node_coordinates[1] = Point{3,Float64}(0.0, 0.0, 0.0)
  node_coordinates[2] = Point{3,Float64}(1.0, 0.0, 0.0)
  node_coordinates[3] = Point{3,Float64}(0.0, 1.0, 0.0)
  node_coordinates[4] = Point{3,Float64}(1.0, 1.0, 0.0)
  node_coordinates[5] = Point{3,Float64}(0.0, 0.0, 1.0)
  node_coordinates[6] = Point{3,Float64}(1.0, 0.0, 1.0)
  node_coordinates[7] = Point{3,Float64}(0.0, 1.0, 1.0)
  node_coordinates[8] = Point{3,Float64}(1.0, 1.0, 1.0)
  node_coordinates[9] = Point{3,Float64}(2.0, 0.0, 0.0)
  node_coordinates[10] = Point{3,Float64}(2.0, 1.0, 0.0)
  node_coordinates[11] = Point{3,Float64}(2.0, 0.0, 1.0)
  node_coordinates[12] = Point{3,Float64}(2.0, 1.0, 1.0)
  polytope = HEX
  scalar_reffe = Gridap.ReferenceFEs.ReferenceFE(polytope, Gridap.ReferenceFEs.lagrangian, Float64, 1)
  cell_types = collect(Fill(1, length(cell_vertex_lids)))
  cell_reffes = [scalar_reffe]
  grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
    cell_vertex_lids,
    cell_reffes,
    cell_types,
    Gridap.Geometry.NonOriented())
  m = Gridap.Geometry.UnstructuredDiscreteModel(grid)
  labels = get_face_labeling(m)
  labels.d_to_dface_to_entity[1] .= 2
  labels.d_to_dface_to_entity[2] .= 2
  labels.d_to_dface_to_entity[3] .= [2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2]
  labels.d_to_dface_to_entity[4] .= 1
  add_tag!(labels, "boundary", [2])
  add_tag!(labels, "interior", [1])
  m
end

with_mpi() do distribute
  ranks = distribute(LinearIndices((MPI.Comm_size(MPI.COMM_WORLD),)))
  coarse_model = init_coarse_model()
  model = OctreeDistributedDiscreteModel(ranks, coarse_model, 1)
end

The error trace for the crash is here:

ERROR: LoadError: AssertionError: sc_array_object.elem_size == sizeof(p8est_tree_t)
Stacktrace:
  [1] p8est_tree_array_index
    @ ~/.julia/packages/P4est_wrapper/tmwT8/src/api_fixes.jl:24 [inlined]
  [2] p8est_tree_array_index
    @ ~/.julia/packages/P4est_wrapper/tmwT8/src/api_fixes.jl:57 [inlined]
  [3] (::GridapP4est.var"#13#14"{3, Ptr{P4est_wrapper.p8est_connectivity}, Vector{Int32}, P4est_wrapper.p8est, P4est_wrapper.p8est_ghost_t, Int64})(cell_vertex_lids::JaggedArray{Int64, Int32}, nl::Int64)
    @ GridapP4est ~/.julia/packages/GridapP4est/MdeY6/src/UniformlyRefinedForestOfOctreesDiscreteModels.jl:376
  [4] map(::GridapP4est.var"#13#14"{3, Ptr{P4est_wrapper.p8est_connectivity}, Vector{Int32}, P4est_wrapper.p8est, P4est_wrapper.p8est_ghost_t, Int64}, ::MPIArray{JaggedArray{Int64, Int32}, 1}, ::MPIArray{Int64, 1})
    @ PartitionedArrays ~/.julia/packages/PartitionedArrays/py6uo/src/mpi_array.jl:229
  [5] generate_node_coordinates(#unused#::Type{Val{3}}, cell_vertex_lids::MPIArray{JaggedArray{Int64, Int32}, 1}, nlvertices::MPIArray{Int64, 1}, ptr_pXest_connectivity::Ptr{P4est_wrapper.p8est_connectivity}, ptr_pXest::Ptr{P4est_wrapper.p8est}, ptr_pXest_ghost::Ptr{P4est_wrapper.p8est_ghost_t})
    @ GridapP4est ~/.julia/packages/GridapP4est/MdeY6/src/UniformlyRefinedForestOfOctreesDiscreteModels.jl:366
  [6] setup_distributed_discrete_model(#unused#::Type{Val{3}}, parts::MPIArray{Int64, 1}, coarse_discrete_model::Gridap.Geometry.UnstructuredDiscreteModel{3, 3, Float64, Gridap.Geometry.NonOriented}, ptr_pXest_connectivity::Ptr{P4est_wrapper.p8est_connectivity}, ptr_pXest::Ptr{P4est_wrapper.p8est}, ptr_pXest_ghost::Ptr{P4est_wrapper.p8est_ghost_t}, ptr_pXest_lnodes::Ptr{P4est_wrapper.p8est_lnodes})
    @ GridapP4est ~/.julia/packages/GridapP4est/MdeY6/src/UniformlyRefinedForestOfOctreesDiscreteModels.jl:927
  [7] OctreeDistributedDiscreteModel(parts::MPIArray{Int64, 1}, coarse_model::Gridap.Geometry.UnstructuredDiscreteModel{3, 3, Float64, Gridap.Geometry.NonOriented}, num_uniform_refinements::Int64)
    @ GridapP4est ~/.julia/packages/GridapP4est/MdeY6/src/OctreeDistributedDiscreteModels.jl:123
  [8] (::var"#3#4")(distribute::PartitionedArrays.var"#88#89"{MPI.Comm, Bool})
    @ Main ~/Downloads/test_p4est_3d_model.jl:49
  [9] with_mpi(f::var"#3#4"; comm::MPI.Comm, duplicate_comm::Bool)
    @ PartitionedArrays ~/.julia/packages/PartitionedArrays/py6uo/src/mpi_array.jl:70
 [10] with_mpi(f::Function)
    @ PartitionedArrays ~/.julia/packages/PartitionedArrays/py6uo/src/mpi_array.jl:64

If the issues are the same, is there any workaround?

[JordiManyer]

Hi @wei3li . What version of p4est are you using? If it is the same error, I have it fixed on a branch of P4est_wrapper.jl. If you want to try it, you can switch to the branch p4est_2.3+

[wei3li]

Hi @JordiManyer, thank you! Branch p4est_2.3+ of P4est_wrapper.jl works for me.