Skip to content

Implement owned cells strategy #15

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Jun 3, 2020
Merged

Implement owned cells strategy #15

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6b2b06c
Implemented OwnedCellsStrategy extension of AssemblyStrategy
amartinhuertas May 27, 2020
ad9319f
Merge branch 'master' of github.com:gridap/GridapDistributed.jl into …
amartinhuertas Jun 1, 2020
259bb4b
+ Addressed @fverdugo remarks for pull request #15:
amartinhuertas Jun 3, 2020
cbe9d3a
Add a function barrier to:
amartinhuertas Jun 3, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
47 changes: 38 additions & 9 deletions src/DistributedAssemblers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,20 +163,19 @@ end
# This typically requires to loop also over ghost cells
struct RowsComputedLocally <: AssemblyStrategy
part::Int
lid_to_gid::Vector{Int}
lid_to_owner::Vector{Int}
gids::IndexSet
end

function Gridap.FESpaces.row_map(a::RowsComputedLocally,row)
a.lid_to_gid[row]
a.gids.lid_to_gid[row]
end

function Gridap.FESpaces.col_map(a::RowsComputedLocally,col)
a.lid_to_gid[col]
a.gids.lid_to_gid[col]
end

function Gridap.FESpaces.row_mask(a::RowsComputedLocally,row)
a.part == a.lid_to_owner[row]
a.part == a.gids.lid_to_owner[row]
end

function Gridap.FESpaces.col_mask(a::RowsComputedLocally,col)
Expand All @@ -186,11 +185,44 @@ end
function RowsComputedLocally(V::DistributedFESpace)
dgids = V.gids
strategies = DistributedData(dgids) do part, gids
RowsComputedLocally(part,gids.lid_to_gid,gids.lid_to_owner)
RowsComputedLocally(part,gids)
end
DistributedAssemblyStrategy(strategies)
end


struct OwnedCellsStrategy <: AssemblyStrategy
part::Int
dof_gids::IndexSet
cell_gids::IndexSet
end

function Gridap.FESpaces.row_map(a::OwnedCellsStrategy,row)
a.dof_gids.lid_to_gid[row]
end

function Gridap.FESpaces.col_map(a::OwnedCellsStrategy,col)
a.dof_gids.lid_to_gid[col]
end

function Gridap.FESpaces.row_mask(a::OwnedCellsStrategy,row)
true
end

function Gridap.FESpaces.col_mask(a::OwnedCellsStrategy,col)
true
end

function OwnedCellsStrategy(M::DistributedDiscreteModel, V::DistributedFESpace)
dcell_gids = M.gids
ddof_gids = V.gids
strategies = DistributedData(ddof_gids,dcell_gids) do part, dof_gids, cell_gids
OwnedCellsStrategy(part,dof_gids,cell_gids)
end
DistributedAssemblyStrategy(strategies)
end


# TODO this assumes that the global matrix type is the same
# as the local one
function Gridap.FESpaces.SparseMatrixAssembler(
Expand All @@ -208,6 +240,3 @@ function Gridap.FESpaces.SparseMatrixAssembler(

DistributedAssembler(matrix_type,vector_type,dtrial,dtest,assems,dstrategy)
end



2 changes: 0 additions & 2 deletions src/DistributedFESpaces.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -250,5 +250,3 @@ function _fill_max_part_around!(lid_to_owner,cell_to_owner,cell_to_lids)
end
end
end


138 changes: 57 additions & 81 deletions src/DistributedTriangulations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,88 +1,64 @@

function remove_ghost_cells(trian::Triangulation,part::Integer,gids::IndexSet)
function filter_cells_when_needed(strategy::AssemblyStrategy, trian::Triangulation)
@abstractmethod
end

function filter_cells_when_needed(strategy::RowsComputedLocally, trian::Triangulation)
trian
end

function filter_cells_when_needed(strategy::OwnedCellsStrategy, trian::Triangulation)
remove_ghost_cells(trian,strategy.part,strategy.cell_gids)
end

function Gridap.Geometry.Triangulation(strategy::AssemblyStrategy,model::DiscreteModel,args...)
trian = Triangulation(model,args...)
filter_cells_when_needed(strategy,trian)
end

function Gridap.Geometry.BoundaryTriangulation(strategy::AssemblyStrategy,model::DiscreteModel,args...)
trian = BoundaryTriangulation(model,args...)
filter_cells_when_needed(strategy,trian)
end

function Gridap.Geometry.SkeletonTriangulation(strategy::AssemblyStrategy,model::DiscreteModel,args...)
trian = SkeletonTriangulation(model,args...)
filter_cells_when_needed(strategy,trian)
end

function remove_ghost_cells(trian::Triangulation, part::Integer, gids::IndexSet)
tcell_to_mcell = get_cell_id(trian)
mcell_to_isowned = gids.lid_to_owner .== part
tcell_to_isowned = reindex(mcell_to_isowned,tcell_to_mcell)
ocell_to_tcell = findall(tcell_to_isowned)
TriangulationPortion(trian,ocell_to_tcell)
ocell_to_tcell =
findall((x) -> (gids.lid_to_owner[x] == part), tcell_to_mcell)
TriangulationPortion(trian, ocell_to_tcell)
end

function include_ghost_cells(trian::TriangulationPortion)
trian.oldtrian
function remove_ghost_cells(
trian::SkeletonTriangulation,
part::Integer,
gids::IndexSet,
)
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 =
_compute_facets_to_old_facets(cell_id_left, cell_id_right, part, gids)
TriangulationPortion(trian, facets_to_old_facets)
end

#
#
#
#struct DistributedTriangulation
# trians::ScatteredVector{<:Triangulation}
#end
#
#function get_distributed_data(dtrian::DistributedTriangulation)
# dtrian.trians
#end
#
#function Gridap.writevtk(dtrian::DistributedTriangulation,filebase::String;cellfields=Dict())
#
# d = Dict(cellfields)
# thevals = values(d)
# thekeys = keys(d)
#
# do_on_parts(dtrian,thevals...) do part, trian, thevals...
# filebase_part = filebase*"_$(part)"
# cellfields = [ k=>v for (k,v) in zip(thekeys, thevals) ]
# writevtk(trian,filebase_part,cellfields=cellfields)
# end
#
#end
function _compute_facets_to_old_facets(cell_id_left, cell_id_right, part, gids)
facets_to_old_facets = eltype(cell_id_right)[]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member Author

@amartinhuertas amartinhuertas Jun 3, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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

for i = 1:length(cell_id_left)
part_left = gids.lid_to_owner[cell_id_left[i]]
part_right = gids.lid_to_owner[cell_id_right[i]]
max_part_id = max(part_left, part_right)
if (max_part_id == part)
push!(facets_to_old_facets, i)
end
end
facets_to_old_facets
end

#
#function Gridap.Triangulation(dmodel::DistributedDiscreteModel,args...)
# comm = get_comm(dmodel)
# trians = ScatteredVector(comm,dmodel.models) do part, model
# Triangulation(model,args...)
# end
# DistributedTriangulation(trians)
#end
#
#function Gridap.BoundaryTriangulation(dmodel::DistributedDiscreteModel,args...)
# comm = get_comm(dmodel)
# trians = ScatteredVector(comm,dmodel.models) do part, model
# BoundaryTriangulation(model,args...)
# end
# DistributedTriangulation(trians)
#end
#
#function Gridap.SkeletonTriangulation(dmodel::DistributedDiscreteModel,args...)
# comm = get_comm(dmodel)
# trians = ScatteredVector(comm,dmodel.models) do part, model
# SkeletonTriangulation(model,args...)
# end
# DistributedTriangulation(trians)
#end
#
#
#function remove_ghost_cells(dtrian::DistributedTriangulation,dmodel)
#
# trians = ScatteredVector(dtrian,dmodel.gids) do part, trian, gids
#
# tcell_to_mcell = get_cell_id(trian)
# mcell_to_isowned = gids.lid_to_owner .== part
# tcell_to_isowned = reindex(mcell_to_isowned,tcell_to_mcell)
# ocell_to_tcell = findall(tcell_to_isowned)
# TriangulationPortion(trian,ocell_to_tcell)
# end
#
# DistributedTriangulation(trians)
#
#end
#
#function include_ghost_cells(dtrian::DistributedTriangulation)
#
# trians = ScatteredVector(dtrian) do part, trian
# trian.oldtrian
# end
#
# DistributedTriangulation(trians)
#end
#
function include_ghost_cells(trian::TriangulationPortion)
trian.oldtrian
end
1 change: 1 addition & 0 deletions src/GridapDistributed.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ export DistributedVector
export exchange!

export RowsComputedLocally
export OwnedCellsStrategy

export remove_ghost_cells
export include_ghost_cells
Expand Down
Loading