WIP: Topologies using CollectionOfVectors

[KnutAM]

Introduces the type CollectionOfVectors which basically allows you to either indexing via a hash like for an OrderedDict or with Ints or CartesianIndex into an Array of predefined size, to get data stored in a linear Vector as a view. Each entry can have different lengths.

Then use this to create topologies, for now demonstrated by the following topologies

• ExclusiveTopology2: Same as ExclusiveTopology, but storage by CollectionOfVectors to reduce the number of small vectors
• CoverTopology: Same as CoverTopology in FerriteDistributed (not tested, but just to check that reuse is possible)
• EntityTopology: Another test where the instances BoundaryIndexes of a unique entity is indexed by looking up the unique entity, e.g. by indexing with the tuple sortface_fast.

(Not all methods implemented yet, but that should be relatively straight forward)

Timing showing stable construction time

grid = generate_grid(Hexahedron, (80, 80, 80));

# ExclusiveTopology2 (pr)
@time Ferrite.ExclusiveTopology2(grid);
# 4.298167 seconds (940.02 k allocations: 2.392 GiB, 10.48% gc time, 15.64% compilation time)
@time Ferrite.ExclusiveTopology2(grid);
# 3.457271 seconds (114 allocations: 2.330 GiB, 7.15% gc time)
@time Ferrite.ExclusiveTopology2(grid);
# 3.385843 seconds (114 allocations: 2.330 GiB, 5.49% gc time)
@time Ferrite.ExclusiveTopology2(grid);
# 3.441457 seconds (114 allocations: 2.330 GiB, 5.87% gc time)

# ExclusiveTopology (master)
@time Ferrite.ExclusiveTopology(grid);
# 8.255457 seconds (33.44 M allocations: 3.703 GiB, 51.47% gc time, 3.15% compilation time)
@time Ferrite.ExclusiveTopology(grid);
# 11.470123 seconds (33.02 M allocations: 3.675 GiB, 66.44% gc time)
@time Ferrite.ExclusiveTopology(grid);
# 23.333075 seconds (33.02 M allocations: 3.675 GiB, 81.44% gc time)

Difference from ArrayOfArrays

The difference from ArrayOfArrays.VectorOfArrays is that the content is always an AbstractVector (specifically a view(::Vector{Float64}, ::UnitRange{Int})), but the outer type can act as a Vector, Matrix, Array, or OrderedDict allowing different indexing.

Old notes

While following the discussion on Slack today, I played around with this idea which doesn't allocate many small vectors. It also uses the unique representations of entities to get the indices.

It is also a bit fast on construction, but slower in usage since it uses hashes to lookup. But the main point is that has much fewer allocations, and ExclusiveTopology brings my computer to a halt due to when creating it twice for a large grid...

julia> grid = generate_grid(Hexahedron, (100,100,100));

julia> @time Ferrite.MaterializedTopology(grid);
  6.347230 seconds (1.31 M allocations: 1.912 GiB, 8.60% gc time, 15.09% compilation time)

julia> @time Ferrite.MaterializedTopology(grid);
  5.112349 seconds (302 allocations: 1.827 GiB, 4.45% gc time)

julia> @time Ferrite.MaterializedTopology(grid);
  5.114099 seconds (302 allocations: 1.827 GiB, 4.57% gc time)

julia> @time ExclusiveTopology(grid);
 17.979864 seconds (65.02 M allocations: 7.223 GiB, 51.87% gc time, 1.43% compilation time)

julia> @time ExclusiveTopology(grid);
214.620891 seconds (64.58 M allocations: 7.195 GiB, 75.27% gc time)

ExclusiveTopology2 and CoverTopology have similar time for construction.

xrefs: #827

[codecov]

Codecov Report

Attention: Patch coverage is 57.09091% with 118 lines in your changes are missing coverage. Please review.

Project coverage is 91.94%. Comparing base (231c8b1) to head (fce4180).
Report is 2 commits behind head on master.

Files	Patch %	Lines
src/Grid/topology2.jl	54.59%	84 Missing ⚠️
src/CollectionOfVectors.jl	61.79%	34 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #965      +/-   ##
==========================================
- Coverage   93.69%   91.94%   -1.75%     
==========================================
  Files          36       38       +2     
  Lines        5438     5713     +275     
==========================================
+ Hits         5095     5253     +158     
- Misses        343      460     +117     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[fredrikekre]

How does this compare to @termi-official s stuff in the distributed repo?

[termi-official]

How does this compare to @termi-official s stuff in the distributed repo?

The topology construction in the distributed implementation (CoverTopology) is very similar to the topology in Ferrite. The standard queries in should be essentially constant time.

[KnutAM]

Edit: Outdated

I did this a bit for fun to see how it would work. The dict-approach is quite nice when building the information, since it makes it quite easy to place everything right.

One idea would be to use MaterializedTopology (perhaps without the compress_data! step) to materialize ExclusiveTopology containing e.g. indices::Matrix{UnitRange{Int}} indexed by (cellid, entitynr). The actual information would then be stored in e.g. data::Vector{BI<:BoundaryIndex}. Since all information is already available the size is known in advance, making construction fast. To query the neighborhood would then just be view(data::Vector{BI}, indices[cellid, entitynr]). This lookup would be faster than MaterializedTopology since it doesn't require hash lookup.

[KnutAM]

As pointed out by @termi-official on Slack, this could also be useful for saving cell states in a continuous vector. Of course, this might also belong in a different package, but that could perhaps be done in the future if we choose to use this strategy.

[koehlerson]

julia> @btime Ferrite.getneighborhood($top_ex,$grid,$(FaceIndex(1,3)))
  2.052 ns (0 allocations: 0 bytes)
1-element Vector{FaceIndex}:
 FaceIndex((2, 5))

julia> @btime Ferrite.getneighborhood($top_ex2,$grid,$(FaceIndex(1,3)))
  2.145 ns (0 allocations: 0 bytes)
1-element view(::Vector{FaceIndex}, 2:2) with eltype FaceIndex:
 FaceIndex((2, 5))

julia> @btime Ferrite.getneighborhood($top_et,$grid,$(FaceIndex(1,3)))
  16.957 ns (0 allocations: 0 bytes)
2-element view(::Vector{FaceIndex}, 3:4) with eltype FaceIndex:
 FaceIndex((1, 3))
 FaceIndex((2, 5))

not sure how much the self mapping in the last case of EntityTopology matters. The old ExclusiveTopology and ExclusiveTopology2 is saving more info like vertex_to_cell and cell_neighbor which naturally needs more allocations. Maybe they don't alter the result too much, but it should be taken into account when comparing the allocation numbers. It seems to me that patching the ExclusiveTopology to ExclusiveTopology2 by using some more appropriate data structure under the hood is favorable, since in the end its the getneighborhood call that is called quite often, but again I'm not sure how much the 10x nanosecond result matters in the end, so maybe it's also worth to compare runs of the topology optimization example or something like that.

[KnutAM]

In any case a fast-lookup datastructure is nice, and any of the proposed constructors are fine IMO since they have a few 100 allocations. Adding the EntityTopology kindof just happened because this was a first idea, but I don't think it is worth changing to that.

Yes, I would also aim to replace ExclusiveTopology by ExclusiveTopology2 if you think that would work well @koehlerson? I like having CollectionOfVectors as a general data-structure and use it here, but would like confirmation that doing that is ok before spending more time on optimizing, unit testing, documentation etc.

I'll make a separate PR for that, just to keep these other "ideas" around.

[KristofferC]

Some random comments / questions:

---

The reason for CollectionOfVectors is to allow for more general outer collections? But couldn't we use ArrayofArrays if we store the general outer collection as a mapping to a linear index that we then use for the vector. So like (in pseudocodeish):

struct CollectionOfVectors
    linear_lookup::Dict{Whatever, Int}
    data::ArrayOfArray{T}
end

Base.getindex(cv::CollectionOfVectors, w::Whatever) = cv.data[cv.linear_lookup[w]]

For matrices, you would only have to store the number of rows to do the linear index computation.

Also, when is a Dict interface needed at all? Are there some sparsity that we didn't exploit before when using a Matrix ?

---

It's a bit hard to get an overview of what is changed here since it isn't a diff for the existing implementation. Is it all new or is most of the previous implementation possible to reuse? Ideally, the data structure should be swappable without too much changes to the implementation of the functions since they should have a similar API?

---

The construction benchmarks looks good but would also be good to see query benchmarks.

[KnutAM]

It's a bit hard to get an overview of what is changed here since it isn't a diff for the existing implementation

Sorry about that: Basically all access function should be almost the same, and most part of the construction algorithms are the same, except that it does add!(x_x_neighbor, value, index) instead of push!(x_x_neibbor[index], value).

The reason for CollectionOfVectors is to allow for more general outer collections? But couldn't we use ArrayofArrays if we store the general outer collection as a mapping to a linear index that we then use for the vector. So like (in pseudocodeish):

For the indexing part, yes that would be similar, but gives an additional redirection AFAIU. But for the construction, ColletionOfVector main feature is that we use the same linear data-vector during construction, which works by allocating a certain length in the vector for a given index (length given by sizehint), and then fill this. If full, it will be moved to the end of the datavector, and twice the space is added, and so it continues. After all data is filled, the data is compressed by filling all empty space by left-shifting the data, such that data is densely stored along our datavector. Then call resize! + sizehint! to free the additional memory.

I could not see how this would be done for ArrayOfArrays, seems to me like you either have to allocate a fixed amount of space for each vector from the beginning, or constructo it from a Vector{Vector}, which defeats the purpose of avoiding all of these many small allocations. But perhaps there is a better way?

Also, when is a Dict interface needed at all? Are there some sparsity that we didn't exploit before when using a Matrix ?

For the topology of different elements, there may be some differences (e.g. Hexahedron with 8 vertices vs Tetrahedron with 4). The initial thought with the EntityTopology was to only store one neighborhood for each unique entity, instead of storing the (same) neighborhood for each BoundaryIndex instance of each entity. However, the lookup costs seems to motivate storing more. (Note that we could store less data by having all items in the index matrix point to the same location in data. For constructing this, the dict would be useful but not required).

There are other use-cases, I think @termi-official mentioned state-variable storage (when you have different cells in the grid) and dof-distribution as places where the Dict-indexing would be beneficial.

So to conclude - if someone has an idea for an efficient construction of VectorOfArrays that is simpler, that would be great.

[KnutAM]

For the accessor time, aren't @koehlerson's timings here what you ask for @KristofferC ?

[koehlerson]

for reproducing the timings:

diff --git a/src/CollectionOfVectors.jl b/src/CollectionOfVectors.jl
index 163cebfa4..93b38a520 100644
--- a/src/CollectionOfVectors.jl
+++ b/src/CollectionOfVectors.jl
@@ -4,6 +4,7 @@ struct CollectionOfVectors{IT<:Union{AbstractArray{UnitRange{Int}}, AbstractDict
 end

 Base.getindex(cv::CollectionOfVectors, index) = view(cv.data, cv.indices[index])
+Base.getindex(cv::CollectionOfVectors, index...) = view(cv.data, cv.indices[index...])
 Base.keys(cv::CollectionOfVectors) = keys(cv.indcies)
 Base.values(cv::CollectionOfVectors) = (view(cv.data, idx) for idx in values(cv.indices))
 nonempty_values(cv::CollectionOfVectors) = (view(cv.data, idx) for idx in values(cv.indices) if !isempty(idx))
diff --git a/src/Grid/topology2.jl b/src/Grid/topology2.jl
index 51118b431..974dbf64d 100644
--- a/src/Grid/topology2.jl
+++ b/src/Grid/topology2.jl
@@ -272,6 +272,14 @@ function ExclusiveTopology2(grid::AbstractGrid{sdim}) where sdim
     return ExclusiveTopology2(vertex_to_cell, cell_neighbor, face_face_neighbor, edge_edge_neighbor, vertex_vertex_neighbor, nothing, nothing, nothing)
 end

+function getneighborhood(top::ExclusiveTopology2, grid::AbstractGrid, faceidx::FaceIndex, include_self=false)
+    if include_self
+        return [top.face_face_neighbor[faceidx[1],faceidx[2]]; faceidx]
+    else
+        return top.face_face_neighbor[faceidx[1],faceidx[2]]
+    end
+end
+

[KnutAM]

Superseded by #974