Implement owned cells strategy #15
Conversation
Codecov Report
@@ Coverage Diff @@
## master #15 +/- ##
==========================================
- Coverage 91.31% 91.02% -0.29%
==========================================
Files 12 12
Lines 426 457 +31
==========================================
+ Hits 389 416 +27
- Misses 37 41 +4
Continue to review full report at Codecov.
|
There was a problem hiding this comment.
@amartinhuertas Thaks for the work!
I have added some comments to be addressed.
src/DistributedFESpaces.jl
Outdated
| vector_type::V | ||
| spaces::DistributedData{<:FESpace} | ||
| gids::DistributedIndexSet | ||
| model::DistributedDiscreteModel |
There was a problem hiding this comment.
Is really necessary to include this field here? In some situations is possible to construct a FE Space without a model (e.g., a DG space from a triangulation). Thus, we need to remove this field.
There was a problem hiding this comment.
No, (I realized after your comment) it is not. I was assuming as a dogma the following statement: "it has to be possible to create any type extension of AssemblyStrategy solely from the local portion of a DistributedFESpace". But this is clearly not the case as the signature of a method that returns a DistributedAssemblyStrategy instance is not fixed, so that we have 100% freedom here. (Note that OwnedCellsStrategy requires access to the corresponding IndexSet portion of the DistributedIndexSet object associated to the distributed model).
src/DistributedTriangulations.jl
Outdated
| function filter_cells_when_needed(strategy::OwnedCellsStrategy, model::DiscreteModel, trian::SkeletonTriangulation) | ||
| topo = get_grid_topology(model) | ||
| D = num_cell_dims(model) | ||
| facets_to_cells = Gridap.Geometry.get_faces(topo,D-1,D) |
There was a problem hiding this comment.
I think this should be implemented as a new method of remove_ghost_cells specialized for SkeletonTriangulation in order to promote code re-usage.
Moreover, perhaps, the implementation can be simplified a bit with something like this (to ckech that it works):
function remove_ghost_cells(trian::SkeletonTriangulation,part,cell_gids)
left = remove_ghost_cells(trian.left,part,cell_gids)
right = reindex(trian.right,left.cell_to_oldcell)
SkeletonTriangulation(left,right)
endThere was a problem hiding this comment.
I think this should be implemented as a new method of remove_ghost_cells specialized for SkeletonTriangulation in order to promote code re-usage.
Done! The old method (interface & rationale) resulted from the fact that I did not realize that I can extract the cells sharing a facet already from SkeletonTriangulation, i.e., I do not need the topology data provided by DiscreteModel.
Moreover, perhaps, the implementation can be simplified a bit with something like this (to ckech that it works):
If I am not missing something, I think the implementation you proposed is not correct (I did not check it). In particular, it leads to an "inconsistent" output SkeletonTriangulation (lets call it output_trian) in the sense that it does NOT fullfill the following assertion: num_cells(output_trian.left) == num_cells(output_trian.right). See also the attached scanned diagram to clarify this point.
CamScanner 06-03-2020 12.31.32.pdf
Do u agree or am I totally confused?
Finally, I am wondering if remove_ghost_cells is the appropiate name for the function at hand. Perhaps something like: build_part_local_integration_triangulation (or something in this line) would be more intention-revealing. What do u think?
There was a problem hiding this comment.
it leads to an "inconsistent" output SkeletonTriangulation (lets call it output_trian) in the sense that it does NOT fullfill the following assertion: num_cells(output_trian.left) == num_cells(output_trian.right).
I would say that this assertion is indeed fulfilled. Anyway, I am not sure if my proposed approach is consistent across subdomains. So, let's take your approach.
There was a problem hiding this comment.
I would say that this assertion is indeed fulfilled.
You are right. It is fulfilled. It is guaranteed by the fact that we re-use left.cell_to_old_cell to create the portion of the "right" boundary triangulation of the original skeleton triangulation.
I missed that code statement, and (wrongly) though that it was actually left = remove_ghost_cells(trian.right,part,cell_gids) (my diagram correspond to this latter scenario).
Anyway, I am not sure if my proposed approach is consistent across subdomains.
I see what you mean. Your approach works as far as for any facet shared among two subdomains, the concept of "left" and "right" cell matches in both subdomains. If it doesnt for a given facet, you end up adding the contribution of that facet twice to the global linear system. If am not wrong, at present, at a given processor, the left cell is the one with smallest LOCAL cell identifier among the two. As, in general, a processor is free to select the local identifier of its cells, this is highly dangerous. Thus, my approach ensures that the problem is solved indepently of how the local cells are numbered.
There was a problem hiding this comment.
BTW, another (minor?) consequence (minor drawback?) of your approach vs mine, is that yours requires left to be a TriangulationPortion, otherwise duck typing fails, while mine does not.
|
@amartinhuertas I have added a second comment. |
* Removed model member variable from DistributedFESpace. * Implemented specialized remove_ghost_cells for SkeletonTriangulation. + Improved performance of the most general version of remove_ghost_cells
Thanks for your comments! I addressed them in commit 259bb4b along with an optimized version of the You can now proceed with the second round of review. |
| cell_id_left = get_cell_id(trian.left) | ||
| cell_id_right = get_cell_id(trian.right) | ||
| @assert length(cell_id_left) == length(cell_id_right) | ||
| facets_to_old_facets = eltype(cell_id_right)[] |
There was a problem hiding this comment.
@amartinhuertas I am not sure if this function is 100% type stable. Perhaps a function barrier for the loop will help. The trian object could be type in-stable (now or in the future) since it is a quite high-level object.
There was a problem hiding this comment.
Damm! I always forget about the issue of type instability. I will definitely add the function barrier. But I do not understand why do u say that trian could be type instable as it is declared to be of type SkeletonTriangulation. trian.left and trian.right types are also known. So I do not see whats wrong in this particular scenario.
There was a problem hiding this comment.
For my understanding, I ran the type inference algorithm (in the debugger), and as I expected, there is no particular type instability (now). Results below. Can you think of a future scenario in which there might be type instability? May be for a different implementation of SkeletonTriangulation ?
debug> InteractiveUtils.@code_warntype remove_ghost_cells(trian,strategy.part,strategy.cell_gids)
Variables
#self#::Core.Compiler.Const(GridapDistributed.remove_ghost_cells, false)
trian::Gridap.Geometry.SkeletonTriangulation{1,2,Gridap.Geometry.GenericBoundaryTriangulation{1,2,Gridap.Geometry.TriangulationPortion{1,2,Gridap.Geometry.UnstructuredGrid{1,2,Float64,false}},Gridap.Geometry.CartesianGrid{2,Float64,typeof(identity)},true,Int64}}
part::Int64
gids::GridapDistributed.IndexSet{Array{Int64,1},Array{Int64,1},Dict{Int64,Int32}}
cell_id_left::Array{Int64,1}
cell_id_right::Array{Int64,1}
facets_to_old_facets::Array{Int64,1}
@_8::Union{Nothing, Tuple{Int64,Int64}}
i::Int64
part_left::Int64
part_right::Int64
max_part_id::Int64
Body::Gridap.Geometry.SkeletonTriangulation{1,2,Gridap.Geometry.TriangulationPortion{1,2,Gridap.Geometry.GenericBoundaryTriangulation{1,2,Gridap.Geometry.TriangulationPortion{1,2,Gridap.Geometry.UnstructuredGrid{1,2,Float64,false}},Gridap.Geometry.CartesianGrid{2,Float64,typeof(identity)},true,Int64}}}
1 ─ Core.NewvarNode(:(facets_to_old_facets))
│ Core.NewvarNode(:(@_8))
│ %3 = Base.getproperty(trian, :left)::Gridap.Geometry.GenericBoundaryTriangulation{1,2,Gridap.Geometry.TriangulationPortion{1,2,Gridap.Geometry.UnstructuredGrid{1,2,Float64,false}},Gridap.Geometry.CartesianGrid{2,Float64,typeof(identity)},true,Int64}
│ (cell_id_left = GridapDistributed.get_cell_id(%3))
│ %5 = Base.getproperty(trian, :right)::Gridap.Geometry.GenericBoundaryTriangulation{1,2,Gridap.Geometry.TriangulationPortion{1,2,Gridap.Geometry.UnstructuredGrid{1,2,Float64,false}},Gridap.Geometry.CartesianGrid{2,Float64,typeof(identity)},true,Int64}
│ (cell_id_right = GridapDistributed.get_cell_id(%5))
│ %7 = GridapDistributed.length(cell_id_left)::Int64
│ %8 = GridapDistributed.length(cell_id_right)::Int64
│ %9 = (%7 == %8)::Bool
└── goto #3 if not %9
2 ─ goto #4
3 ─ %12 = Base.AssertionError("length(cell_id_left) == length(cell_id_right)")::AssertionError
└── Base.throw(%12)
4 ┄ %14 = GridapDistributed.eltype(cell_id_right)::Core.Compiler.Const(Int64, false)
│ (facets_to_old_facets = Base.getindex(%14))
│ %16 = GridapDistributed.length(cell_id_left)::Int64
│ %17 = (1:%16)::Core.Compiler.PartialStruct(UnitRange{Int64}, Any[Core.Compiler.Const(1, false), Int64])
│ (@_8 = Base.iterate(%17))
│ %19 = (@_8 === nothing)::Bool
│ %20 = Base.not_int(%19)::Bool
└── goto #9 if not %20
5 ┄ %22 = @_8::Tuple{Int64,Int64}::Tuple{Int64,Int64}
│ (i = Core.getfield(%22, 1))
│ %24 = Core.getfield(%22, 2)::Int64
│ %25 = Base.getproperty(gids, :lid_to_owner)::Array{Int64,1}
│ %26 = Base.getindex(cell_id_left, i)::Int64
│ (part_left = Base.getindex(%25, %26))
│ %28 = Base.getproperty(gids, :lid_to_owner)::Array{Int64,1}
│ %29 = Base.getindex(cell_id_right, i)::Int64
│ (part_right = Base.getindex(%28, %29))
│ (max_part_id = GridapDistributed.max(part_left, part_right))
│ %32 = (max_part_id == part)::Bool
└── goto #7 if not %32
6 ─ GridapDistributed.push!(facets_to_old_facets, i)
7 ┄ (@_8 = Base.iterate(%17, %24))
│ %36 = (@_8 === nothing)::Bool
│ %37 = Base.not_int(%36)::Bool
└── goto #9 if not %37
8 ─ goto #5
9 ┄ %40 = GridapDistributed.TriangulationPortion(trian, facets_to_old_facets)::Gridap.Geometry.SkeletonTriangulation{1,2,Gridap.Geometry.TriangulationPortion{1,2,Gridap.Geometry.GenericBoundaryTriangulation{1,2,Gridap.Geometry.TriangulationPortion{1,2,Gridap.Geometry.UnstructuredGrid{1,2,Float64,false}},Gridap.Geometry.CartesianGrid{2,Float64,typeof(identity)},true,Int64}}}
└── return %40
There was a problem hiding this comment.
trian.left or trian.right can potentially be a type in-stable specialization of Triangulation. For now it is not the case, but who knows in the future...
There was a problem hiding this comment.
trian.left or trian.right can potentially be a type in-stable specialization of Triangulation. For now it is not the case, but who knows in the future...
Ok. This is what I was thinking. Thanks for confirming!
There was a problem hiding this comment.
@amartinhuertas I am not sure if this function is 100% type stable. Perhaps a function barrier for the loop will help. The trian object could be type in-stable (now or in the future) since it is a quite high-level object.
Done in cbe9d3a
function remove_ghost_cells(
trian::SkeletonTriangulation,
part::Integer,
gids::IndexSet,
)
with type stability in mind
I already addressed it. You can proceed with the second round of review. Thanks for your insightful comments that helped a lot to improve the code quality! Finally, any wording from yout side on this remark that I did before?
|
The current name is shorter and already quite reveling. |
No description provided.