From 6b2b06c04920b0f741d9478450086e209cb5d212 Mon Sep 17 00:00:00 2001 From: amartin Date: Wed, 27 May 2020 22:30:25 +1000 Subject: [PATCH 1/3] Implemented OwnedCellsStrategy extension of AssemblyStrategy --- src/DistributedAssemblers.jl | 47 +++++-- src/DistributedFESpaces.jl | 5 +- src/DistributedTriangulations.jl | 124 +++++++----------- src/GridapDistributed.jl | 1 + test/DistributedPoissonDGTests.jl | 201 +++++++++++++++--------------- test/DistributedPoissonTests.jl | 149 +++++++++++----------- 6 files changed, 268 insertions(+), 259 deletions(-) diff --git a/src/DistributedAssemblers.jl b/src/DistributedAssemblers.jl index 02ecbeae..c429b480 100644 --- a/src/DistributedAssemblers.jl +++ b/src/DistributedAssemblers.jl @@ -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) @@ -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(V::DistributedFESpace) + ddof_gids = V.gids + dcell_gids = V.model.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( @@ -208,6 +240,3 @@ function Gridap.FESpaces.SparseMatrixAssembler( DistributedAssembler(matrix_type,vector_type,dtrial,dtest,assems,dstrategy) end - - - diff --git a/src/DistributedFESpaces.jl b/src/DistributedFESpaces.jl index 1cf24a63..1cef4fdc 100644 --- a/src/DistributedFESpaces.jl +++ b/src/DistributedFESpaces.jl @@ -2,6 +2,7 @@ struct DistributedFESpace{V} <: FESpace vector_type::V spaces::DistributedData{<:FESpace} gids::DistributedIndexSet + model::DistributedDiscreteModel end function get_distributed_data(dspace::DistributedFESpace) @@ -80,7 +81,7 @@ function Gridap.TrialFESpace(V::DistributedFESpace,args...) spaces = DistributedData(V.spaces) do part, space TrialFESpace(space,args...) end - DistributedFESpace(V.vector_type,spaces,V.gids) + DistributedFESpace(V.vector_type,spaces,V.gids,V.model) end function Gridap.FESpace(::Type{V};model::DistributedDiscreteModel,kwargs...) where V @@ -180,7 +181,7 @@ function DistributedFESpace(::Type{V}; model::DistributedDiscreteModel,kwargs... gids = DistributedIndexSet(init_free_gids,comm,ngids, part_to_lid_to_gid,part_to_lid_to_owner,part_to_ngids) - DistributedFESpace(V,spaces,gids) + DistributedFESpace(V,spaces,gids,model) end function _update_lid_to_gid!(lid_to_gid,cell_to_lids,cell_to_gids,cell_to_owner,lid_to_owner) diff --git a/src/DistributedTriangulations.jl b/src/DistributedTriangulations.jl index 4ea57e2e..0a96952a 100644 --- a/src/DistributedTriangulations.jl +++ b/src/DistributedTriangulations.jl @@ -1,4 +1,52 @@ +function filter_cells_when_needed(strategy::AssemblyStrategy, model::DiscreteModel, trian::Triangulation) + @abstractmethod +end + +function filter_cells_when_needed(strategy::RowsComputedLocally, model::DiscreteModel, trian::Triangulation) + trian +end + +function filter_cells_when_needed(strategy::OwnedCellsStrategy, model::DiscreteModel, trian::Triangulation) + remove_ghost_cells(trian,strategy.part,strategy.cell_gids) +end + +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) + facet_lids = Gridap.Geometry.get_cell_id(trian.left.face_trian) + facets_to_old_facets = eltype(facet_lids)[] + for (i,facet_lid) in enumerate(facet_lids) + max_part_id=-1 + for j=facets_to_cells.ptrs[facet_lid]:facets_to_cells.ptrs[facet_lid+1]-1 + cell_lid = facets_to_cells.data[j] + cell_part = strategy.cell_gids.lid_to_owner[cell_lid] + max_part_id=max(max_part_id, cell_part) + end + if (strategy.part == max_part_id) + push!(facets_to_old_facets, i) + end + end + TriangulationPortion(trian,facets_to_old_facets) + #trian +end + +function Gridap.Geometry.Triangulation(strategy::AssemblyStrategy,model::DiscreteModel,args...) + trian = Triangulation(model,args...) + filter_cells_when_needed(strategy,model,trian) +end + +function Gridap.Geometry.BoundaryTriangulation(strategy::AssemblyStrategy,model::DiscreteModel,args...) + trian = BoundaryTriangulation(model,args...) + filter_cells_when_needed(strategy,model,trian) +end + +function Gridap.Geometry.SkeletonTriangulation(strategy::AssemblyStrategy,model::DiscreteModel,args...) + trian = SkeletonTriangulation(model,args...) + filter_cells_when_needed(strategy,model,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 @@ -10,79 +58,3 @@ end function include_ghost_cells(trian::TriangulationPortion) trian.oldtrian 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 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 -# diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl index 8c0f2c75..127f4628 100644 --- a/src/GridapDistributed.jl +++ b/src/GridapDistributed.jl @@ -32,6 +32,7 @@ export DistributedVector export exchange! export RowsComputedLocally +export OwnedCellsStrategy export remove_ghost_cells export include_ghost_cells diff --git a/test/DistributedPoissonDGTests.jl b/test/DistributedPoissonDGTests.jl index 11764da3..19dbafc7 100644 --- a/test/DistributedPoissonDGTests.jl +++ b/test/DistributedPoissonDGTests.jl @@ -6,105 +6,108 @@ using Gridap.FESpaces using GridapDistributed using SparseArrays -# Select matrix and vector types for discrete problem -# Note that here we use serial vectors and matrices -# but the assembly is distributed -T = Float64 -vector_type = Vector{T} -matrix_type = SparseMatrixCSC{T,Int} - -# Manufactured solution -u(x) = x[1]*(x[1]-1)*x[2]*(x[2]-1) -f(x) = - Δ(u)(x) -ud(x) = zero(x[1]) - -# Discretization -subdomains = (2,2) -domain = (0,1,0,1) -cells = (4,4) -comm = SequentialCommunicator(subdomains) -model = CartesianDiscreteModel(comm,subdomains,domain,cells) -const h = (domain[2]-domain[1]) / cells[1] - -# FE Spaces -order = 2 -V = FESpace( - vector_type, valuetype=Float64, reffe=:Lagrangian, order=order, - model=model, conformity=:L2) - -U = TrialFESpace(V) - -# Terms in the weak form -terms = DistributedData(model) do part, (model,gids) - - γ = 10 - degree = 2*order - - trian = get_triangulation(model) - btrian = BoundaryTriangulation(model) - strian = SkeletonTriangulation(model) - - quad = CellQuadrature(trian,degree) - bquad = CellQuadrature(btrian,degree) - squad = CellQuadrature(strian,degree) - - bn = get_normal_vector(btrian) - sn = get_normal_vector(strian) - - a(u,v) = inner(∇(v),∇(u)) - l(v) = v*f - t_Ω = AffineFETerm(a,l,trian,quad) - - a_Γd(u,v) = (γ/h)*v*u - v*(bn*∇(u)) - (bn*∇(v))*u - l_Γd(v) = (γ/h)*v*ud - (bn*∇(v))*ud - t_Γd = AffineFETerm(a_Γd,l_Γd,btrian,bquad) - - a_Γ(u,v) = (γ/h)*jump(v*sn)*jump(u*sn) - jump(v*sn)*mean(∇(u)) - mean(∇(v))*jump(u*sn) - t_Γ = LinearFETerm(a_Γ,strian,squad) - - (t_Ω,t_Γ,t_Γd) +function run(assembly_strategy::AbstractString) + # Select matrix and vector types for discrete problem + # Note that here we use serial vectors and matrices + # but the assembly is distributed + T = Float64 + vector_type = Vector{T} + matrix_type = SparseMatrixCSC{T,Int} + + # Manufactured solution + u(x) = x[1] * (x[1] - 1) * x[2] * (x[2] - 1) + f(x) = -Δ(u)(x) + ud(x) = zero(x[1]) + + # Discretization + subdomains = (2, 2) + domain = (0, 1, 0, 1) + cells = (4, 4) + comm = SequentialCommunicator(subdomains) + model = CartesianDiscreteModel(comm, subdomains, domain, cells) + h = (domain[2] - domain[1]) / cells[1] + + # FE Spaces + order = 2 + V = FESpace( + vector_type, + valuetype = Float64, + reffe = :Lagrangian, + order = order, + model = model, + conformity = :L2, + ) + + U = TrialFESpace(V) + + if (assembly_strategy == "RowsComputedLocally") + strategy = RowsComputedLocally(V) + elseif (assembly_strategy == "OwnedCellsStrategy") + strategy = OwnedCellsStrategy(V) + else + @assert false "Unknown AssemblyStrategy: $(assembly_strategy)" + end + + # Terms in the weak form + terms = DistributedData(model, strategy) do part, (model, gids), strategy + γ = 10 + degree = 2 * order + + trian = Triangulation(strategy, model) + btrian = BoundaryTriangulation(strategy, model) + strian = SkeletonTriangulation(strategy, model) + + quad = CellQuadrature(trian, degree) + bquad = CellQuadrature(btrian, degree) + squad = CellQuadrature(strian, degree) + + bn = get_normal_vector(btrian) + sn = get_normal_vector(strian) + + a(u, v) = inner(∇(v), ∇(u)) + l(v) = v * f + t_Ω = AffineFETerm(a, l, trian, quad) + + a_Γd(u, v) = (γ / h) * v * u - v * (bn * ∇(u)) - (bn * ∇(v)) * u + l_Γd(v) = (γ / h) * v * ud - (bn * ∇(v)) * ud + t_Γd = AffineFETerm(a_Γd, l_Γd, btrian, bquad) + + a_Γ(u, v) = (γ / h) * jump(v * sn) * jump(u * sn) - jump(v * sn) * mean(∇(u)) - + mean(∇(v)) * jump(u * sn) + t_Γ = LinearFETerm(a_Γ, strian, squad) + (t_Ω, t_Γ, t_Γd) + end + + # Assembly + assem = SparseMatrixAssembler(matrix_type, vector_type, U, V, strategy) + op = AffineFEOperator(assem, terms) + + # FE solution + uh = solve(op) + + # Error norms and print solution + sums = DistributedData(model, uh) do part, (model, gids), uh + trian = Triangulation(model) + owned_trian = remove_ghost_cells(trian, part, gids) + owned_quad = CellQuadrature(owned_trian, 2 * order) + owned_uh = restrict(uh, owned_trian) + writevtk(owned_trian, "results_$part", cellfields = ["uh" => owned_uh]) + e = u - owned_uh + l2(u) = u * u + h1(u) = u * u + ∇(u) * ∇(u) + e_l2 = sum(integrate(l2(e), owned_trian, owned_quad)) + e_h1 = sum(integrate(h1(e), owned_trian, owned_quad)) + e_l2, e_h1 + end + e_l2_h1 = gather(sums) + e_l2 = sum(map(i -> i[1], e_l2_h1)) + e_h1 = sum(map(i -> i[2], e_l2_h1)) + tol = 1.0e-9 + @test e_l2 < tol + @test e_h1 < tol end -# Chose parallel assembly strategy -strategy = RowsComputedLocally(V) +run("RowsComputedLocally") +run("OwnedCellsStrategy") -# Assembly -assem = SparseMatrixAssembler(matrix_type, vector_type, U, V, strategy) -op = AffineFEOperator(assem,terms) - -# FE solution -op = AffineFEOperator(assem,terms) -uh = solve(op) - -# Error norms and print solution -sums = DistributedData(model,uh) do part, (model,gids), uh - - trian = Triangulation(model) - - owned_trian = remove_ghost_cells(trian,part,gids) - owned_quad = CellQuadrature(owned_trian,2*order) - owned_uh = restrict(uh,owned_trian) - - writevtk(owned_trian,"results_$part",cellfields=["uh"=>owned_uh]) - - e = u - owned_uh - - l2(u) = u*u - h1(u) = u*u + ∇(u)*∇(u) - - e_l2 = sum(integrate(l2(e),owned_trian,owned_quad)) - e_h1 = sum(integrate(h1(e),owned_trian,owned_quad)) - - e_l2, e_h1 -end - -e_l2_h1 = gather(sums) - -e_l2 = sum(map(i->i[1],e_l2_h1)) -e_h1 = sum(map(i->i[2],e_l2_h1)) - -tol = 1.0e-9 -@test e_l2 < tol -@test e_h1 < tol - -end # module +end # module# diff --git a/test/DistributedPoissonTests.jl b/test/DistributedPoissonTests.jl index fa68d350..b184ab16 100644 --- a/test/DistributedPoissonTests.jl +++ b/test/DistributedPoissonTests.jl @@ -6,80 +6,83 @@ using Gridap.FESpaces using GridapDistributed using SparseArrays -# Select matrix and vector types for discrete problem -# Note that here we use serial vectors and matrices -# but the assembly is distributed -T = Float64 -vector_type = Vector{T} -matrix_type = SparseMatrixCSC{T,Int} - -# Manufactured solution -u(x) = x[1] + x[2] -f(x) = - Δ(u)(x) - -# Discretization -subdomains = (2,2) -domain = (0,1,0,1) -cells = (4,4) -comm = SequentialCommunicator(subdomains) -model = CartesianDiscreteModel(comm,subdomains,domain,cells) - -# FE Spaces -order = 1 -V = FESpace( - vector_type, valuetype=Float64, reffe=:Lagrangian, order=order, - model=model, conformity=:H1, dirichlet_tags="boundary") - -U = TrialFESpace(V,u) - -# Terms in the weak form -terms = DistributedData(model) do part, (model,gids) - - trian = Triangulation(model) - - degree = 2*order - quad = CellQuadrature(trian,degree) - - a(u,v) = ∇(v)*∇(u) - l(v) = v*f - t1 = AffineFETerm(a,l,trian,quad) - - (t1,) -end - -# Chose parallel assembly strategy -strategy = RowsComputedLocally(V) - -# Assembler -assem = SparseMatrixAssembler(matrix_type, vector_type, U, V, strategy) - -# FE solution -op = AffineFEOperator(assem,terms) -uh = solve(op) - -# Error norms and print solution -sums = DistributedData(model,uh) do part, (model,gids), uh - - trian = Triangulation(model) - - owned_trian = remove_ghost_cells(trian,part,gids) - owned_quad = CellQuadrature(owned_trian,2*order) - owned_uh = restrict(uh,owned_trian) - - writevtk(owned_trian,"results_$part",cellfields=["uh"=>owned_uh]) - - e = u - owned_uh - - l2(u) = u*u - - sum(integrate(l2(e),owned_trian,owned_quad)) - +function run(assembly_strategy::AbstractString) + # Select matrix and vector types for discrete problem + # Note that here we use serial vectors and matrices + # but the assembly is distributed + T = Float64 + vector_type = Vector{T} + matrix_type = SparseMatrixCSC{T,Int} + + # Manufactured solution + u(x) = x[1] + x[2] + f(x) = -Δ(u)(x) + + # Discretization + subdomains = (2, 2) + domain = (0, 1, 0, 1) + cells = (4, 4) + comm = SequentialCommunicator(subdomains) + model = CartesianDiscreteModel(comm, subdomains, domain, cells) + + # FE Spaces + order = 1 + V = FESpace( + vector_type, + valuetype = Float64, + reffe = :Lagrangian, + order = order, + model = model, + conformity = :H1, + dirichlet_tags = "boundary", + ) + + U = TrialFESpace(V, u) + + if (assembly_strategy == "RowsComputedLocally") + strategy = RowsComputedLocally(V) + elseif (assembly_strategy == "OwnedCellsStrategy") + strategy = OwnedCellsStrategy(V) + else + @assert false "Unknown AssemblyStrategy: $(assembly_strategy)" + end + + # Terms in the weak form + terms = DistributedData(model, strategy) do part, (model, gids), strategy + trian = Triangulation(strategy, model) + degree = 2 * order + quad = CellQuadrature(trian, degree) + a(u, v) = ∇(v) * ∇(u) + l(v) = v * f + t1 = AffineFETerm(a, l, trian, quad) + (t1,) + end + + # Assembler + assem = SparseMatrixAssembler(matrix_type, vector_type, U, V, strategy) + + # FE solution + op = AffineFEOperator(assem, terms) + uh = solve(op) + + # Error norms and print solution + sums = DistributedData(model, uh) do part, (model, gids), uh + trian = Triangulation(model) + owned_trian = remove_ghost_cells(trian, part, gids) + owned_quad = CellQuadrature(owned_trian, 2 * order) + owned_uh = restrict(uh, owned_trian) + writevtk(owned_trian, "results_$part", cellfields = ["uh" => owned_uh]) + e = u - owned_uh + l2(u) = u * u + sum(integrate(l2(e), owned_trian, owned_quad)) + end + e_l2 = sum(gather(sums)) + + tol = 1.0e-9 + @test e_l2 < tol end -e_l2 = sum(gather(sums)) - -tol = 1.0e-9 -@test e_l2 < tol - +run("RowsComputedLocally") +run("OwnedCellsStrategy") end # module From 259bb4b2e7d7a62f89386362b19b44d8c7dc7a8b Mon Sep 17 00:00:00 2001 From: amartin Date: Wed, 3 Jun 2020 13:28:13 +1000 Subject: [PATCH 2/3] + Addressed @fverdugo remarks for pull request #15: * Removed model member variable from DistributedFESpace. * Implemented specialized remove_ghost_cells for SkeletonTriangulation. + Improved performance of the most general version of remove_ghost_cells --- src/DistributedAssemblers.jl | 4 +- src/DistributedFESpaces.jl | 7 +--- src/DistributedTriangulations.jl | 62 +++++++++++++++---------------- test/DistributedPoissonDGTests.jl | 2 +- test/DistributedPoissonTests.jl | 2 +- 5 files changed, 36 insertions(+), 41 deletions(-) diff --git a/src/DistributedAssemblers.jl b/src/DistributedAssemblers.jl index c429b480..fde7f601 100644 --- a/src/DistributedAssemblers.jl +++ b/src/DistributedAssemblers.jl @@ -213,9 +213,9 @@ function Gridap.FESpaces.col_mask(a::OwnedCellsStrategy,col) true end -function OwnedCellsStrategy(V::DistributedFESpace) +function OwnedCellsStrategy(M::DistributedDiscreteModel, V::DistributedFESpace) + dcell_gids = M.gids ddof_gids = V.gids - dcell_gids = V.model.gids strategies = DistributedData(ddof_gids,dcell_gids) do part, dof_gids, cell_gids OwnedCellsStrategy(part,dof_gids,cell_gids) end diff --git a/src/DistributedFESpaces.jl b/src/DistributedFESpaces.jl index 1cef4fdc..ea9c7150 100644 --- a/src/DistributedFESpaces.jl +++ b/src/DistributedFESpaces.jl @@ -2,7 +2,6 @@ struct DistributedFESpace{V} <: FESpace vector_type::V spaces::DistributedData{<:FESpace} gids::DistributedIndexSet - model::DistributedDiscreteModel end function get_distributed_data(dspace::DistributedFESpace) @@ -81,7 +80,7 @@ function Gridap.TrialFESpace(V::DistributedFESpace,args...) spaces = DistributedData(V.spaces) do part, space TrialFESpace(space,args...) end - DistributedFESpace(V.vector_type,spaces,V.gids,V.model) + DistributedFESpace(V.vector_type,spaces,V.gids) end function Gridap.FESpace(::Type{V};model::DistributedDiscreteModel,kwargs...) where V @@ -181,7 +180,7 @@ function DistributedFESpace(::Type{V}; model::DistributedDiscreteModel,kwargs... gids = DistributedIndexSet(init_free_gids,comm,ngids, part_to_lid_to_gid,part_to_lid_to_owner,part_to_ngids) - DistributedFESpace(V,spaces,gids,model) + DistributedFESpace(V,spaces,gids) end function _update_lid_to_gid!(lid_to_gid,cell_to_lids,cell_to_gids,cell_to_owner,lid_to_owner) @@ -251,5 +250,3 @@ function _fill_max_part_around!(lid_to_owner,cell_to_owner,cell_to_lids) end end end - - diff --git a/src/DistributedTriangulations.jl b/src/DistributedTriangulations.jl index 0a96952a..d2ca979e 100644 --- a/src/DistributedTriangulations.jl +++ b/src/DistributedTriangulations.jl @@ -1,58 +1,56 @@ -function filter_cells_when_needed(strategy::AssemblyStrategy, model::DiscreteModel, trian::Triangulation) +function filter_cells_when_needed(strategy::AssemblyStrategy, trian::Triangulation) @abstractmethod end -function filter_cells_when_needed(strategy::RowsComputedLocally, model::DiscreteModel, trian::Triangulation) +function filter_cells_when_needed(strategy::RowsComputedLocally, trian::Triangulation) trian end -function filter_cells_when_needed(strategy::OwnedCellsStrategy, model::DiscreteModel, trian::Triangulation) +function filter_cells_when_needed(strategy::OwnedCellsStrategy, trian::Triangulation) remove_ghost_cells(trian,strategy.part,strategy.cell_gids) end -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) - facet_lids = Gridap.Geometry.get_cell_id(trian.left.face_trian) - facets_to_old_facets = eltype(facet_lids)[] - for (i,facet_lid) in enumerate(facet_lids) - max_part_id=-1 - for j=facets_to_cells.ptrs[facet_lid]:facets_to_cells.ptrs[facet_lid+1]-1 - cell_lid = facets_to_cells.data[j] - cell_part = strategy.cell_gids.lid_to_owner[cell_lid] - max_part_id=max(max_part_id, cell_part) - end - if (strategy.part == max_part_id) - push!(facets_to_old_facets, i) - end - end - TriangulationPortion(trian,facets_to_old_facets) - #trian -end - function Gridap.Geometry.Triangulation(strategy::AssemblyStrategy,model::DiscreteModel,args...) trian = Triangulation(model,args...) - filter_cells_when_needed(strategy,model,trian) + 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,model,trian) + 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,model,trian) + filter_cells_when_needed(strategy,trian) end -function remove_ghost_cells(trian::Triangulation,part::Integer,gids::IndexSet) +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 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 = eltype(cell_id_right)[] + 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 + TriangulationPortion(trian, facets_to_old_facets) end function include_ghost_cells(trian::TriangulationPortion) diff --git a/test/DistributedPoissonDGTests.jl b/test/DistributedPoissonDGTests.jl index 19dbafc7..e21a03f4 100644 --- a/test/DistributedPoissonDGTests.jl +++ b/test/DistributedPoissonDGTests.jl @@ -43,7 +43,7 @@ function run(assembly_strategy::AbstractString) if (assembly_strategy == "RowsComputedLocally") strategy = RowsComputedLocally(V) elseif (assembly_strategy == "OwnedCellsStrategy") - strategy = OwnedCellsStrategy(V) + strategy = OwnedCellsStrategy(model,V) else @assert false "Unknown AssemblyStrategy: $(assembly_strategy)" end diff --git a/test/DistributedPoissonTests.jl b/test/DistributedPoissonTests.jl index b184ab16..3a097e4f 100644 --- a/test/DistributedPoissonTests.jl +++ b/test/DistributedPoissonTests.jl @@ -42,7 +42,7 @@ function run(assembly_strategy::AbstractString) if (assembly_strategy == "RowsComputedLocally") strategy = RowsComputedLocally(V) elseif (assembly_strategy == "OwnedCellsStrategy") - strategy = OwnedCellsStrategy(V) + strategy = OwnedCellsStrategy(model,V) else @assert false "Unknown AssemblyStrategy: $(assembly_strategy)" end From cbe9d3a22fe3937b95dbe10f4331e357c8db6355 Mon Sep 17 00:00:00 2001 From: amartin Date: Wed, 3 Jun 2020 18:16:49 +1000 Subject: [PATCH 3/3] Add a function barrier to: function remove_ghost_cells( trian::SkeletonTriangulation, part::Integer, gids::IndexSet, ) with type stability in mind --- src/DistributedTriangulations.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/DistributedTriangulations.jl b/src/DistributedTriangulations.jl index d2ca979e..7d111832 100644 --- a/src/DistributedTriangulations.jl +++ b/src/DistributedTriangulations.jl @@ -41,6 +41,12 @@ function remove_ghost_cells( 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 + +function _compute_facets_to_old_facets(cell_id_left, cell_id_right, part, gids) facets_to_old_facets = eltype(cell_id_right)[] for i = 1:length(cell_id_left) part_left = gids.lid_to_owner[cell_id_left[i]] @@ -50,7 +56,7 @@ function remove_ghost_cells( push!(facets_to_old_facets, i) end end - TriangulationPortion(trian, facets_to_old_facets) + facets_to_old_facets end function include_ghost_cells(trian::TriangulationPortion)