Parallelization

Project.toml

@@ -8,12 +8,16 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
+GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+GridapDistributed = "0.2.6"
 DrWatson = "1"
 Gridap = "0.17"

src/Distributed/HybridAffineFEOperators.jl

@@ -0,0 +1,85 @@
+function _merge_bulk_and_skeleton_contributions(
+    matcontribs::GridapDistributed.DistributedDomainContribution,
+    veccontribs::GridapDistributed.DistributedDomainContribution)
+  dmat,dvec=map_parts(matcontribs.contribs,veccontribs.contribs) do matseq, vecseq
+    _merge_bulk_and_skeleton_contributions(matseq,vecseq)
+  end
+  GridapDistributed.DistributedDomainContribution(dmat),
+          GridapDistributed.DistributedDomainContribution(dvec)
+end
+
+function Gridap.FESpaces._pair_contribution_when_possible(
+    obiform::GridapDistributed.DistributedDomainContribution,
+    oliform::GridapDistributed.DistributedDomainContribution)
+  dmv,dm,dv=map_parts(obiform.contribs,oliform.contribs) do obiform, oliform
+    Gridap.FESpaces._pair_contribution_when_possible(obiform,oliform)
+  end
+  GridapDistributed.DistributedDomainContribution(dmv),
+    GridapDistributed.DistributedDomainContribution(dm),
+      GridapDistributed.DistributedDomainContribution(dv)
+end
+
+function Gridap.FESpaces._collect_cell_matrix_and_vector(
+  trial  :: GridapDistributed.DistributedFESpace,
+  test   :: GridapDistributed.DistributedFESpace,
+  matvec :: GridapDistributed.DistributedDomainContribution,
+  mat    :: GridapDistributed.DistributedDomainContribution,
+  vec    :: GridapDistributed.DistributedDomainContribution)
+  map_parts(GridapDistributed.local_views(trial),
+            GridapDistributed.local_views(test),
+            matvec.contribs,mat.contribs,vec.contribs) do trial, test, matvec, mat, vec
+     Gridap.FESpaces._collect_cell_matrix_and_vector(trial,test,matvec,mat,vec)
+  end
+end
+
+function _add_static_condensation(
+   matvec::GridapDistributed.DistributedDomainContribution,
+   bulk_fields, skeleton_fields)
+   dmv=map_parts(matvec.contribs) do matvec
+    _add_static_condensation(matvec,bulk_fields,skeleton_fields)
+  end
+  GridapDistributed.DistributedDomainContribution(dmv)
+end
+
+function Gridap.FESpaces._attach_dirichlet(
+  matvec:: GridapDistributed.DistributedDomainContribution,
+  mat   :: GridapDistributed.DistributedDomainContribution,
+  uhd)
+  dmv,dm=map_parts(matvec.contribs,mat.contribs,uhd.fields) do matvec, mat, uhd
+    m,v=Gridap.FESpaces._attach_dirichlet(matvec,mat, uhd)
+    m,v
+  end
+  GridapDistributed.DistributedDomainContribution(dmv),
+      GridapDistributed.DistributedDomainContribution(dm)
+end
+
+
+function _compute_hybridizable_from_skeleton_free_dof_values(
+    skeleton_fe_function :: GridapDistributed.DistributedCellField,
+    trial_hybridizable   :: GridapDistributed.DistributedFESpace,
+    test_hybridizable    :: GridapDistributed.DistributedFESpace,
+    trial_skeleton       :: GridapDistributed.DistributedFESpace,
+    matvec               :: GridapDistributed.DistributedDomainContribution,
+    bulk_fields, skeleton_fields)
+
+    function f(skel_fe, trial_hyb, test_hyb, trial_skel, matvec)
+      _compute_hybridizable_from_skeleton_free_dof_values(skel_fe,
+                                                          trial_hyb,
+                                                          test_hyb,
+                                                          trial_skel,
+                                                          matvec,
+                                                          bulk_fields,
+                                                          skeleton_fields)
+    end
+    # f(skeleton_fe_function.fields.parts[1],
+    #   trial_hybridizable.part_fe_space.parts[1],
+    #   test_hybridizable.part_fe_space.parts[1],
+    #   trial_skeleton.spaces.parts[1],
+    #   matvec.contribs.parts[1])
+    values=map_parts(f,GridapDistributed.local_views(skeleton_fe_function),
+             GridapDistributed.local_views(trial_hybridizable),
+             GridapDistributed.local_views(test_hybridizable),
+             GridapDistributed.local_views(trial_skeleton),
+             GridapDistributed.local_views(matvec))
+    PVector(values,trial_hybridizable.gids)
+end

src/Distributed/Skeleton.jl

@@ -0,0 +1,49 @@
+
+function SkeletonTriangulation(model::GridapDistributed.DistributedDiscreteModel)
+  SkeletonTriangulation(no_ghost,model::GridapDistributed.DistributedDiscreteModel)
+end
+
+
+function _sign_flips(model::GridapDistributed.DistributedDiscreteModel)
+  cell_reffes = map_parts(model.models) do m
+    basis,reffe_args,reffe_kwargs = ReferenceFE(raviart_thomas,Float64,0)
+    reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...)
+  end
+  GridapDistributed._generate_sign_flips(model,cell_reffes)
+end
+
+function SkeletonTriangulation(portion::GridapDistributed.WithGhost,
+                               model::GridapDistributed.DistributedDiscreteModel)
+   dtrians=map_parts(model.models,_sign_flips(model)) do model, sign_flip
+      Skeleton(model,sign_flip)
+   end
+   GridapDistributed.DistributedTriangulation(dtrians,model)
+end
+
+function SkeletonTriangulation(portion::GridapDistributed.NoGhost,
+                               model::GridapDistributed.DistributedDiscreteModel)
+  cell_gids=get_cell_gids(model)
+  dtrians=map_parts(model.models,
+                    _sign_flips(model),
+                    cell_gids.partition) do model, sign_flip, partition
+    cell_to_parent_cell=
+      findall([partition.lid_to_part[cell]==partition.part
+                 for cell=1:length(partition.lid_to_part)])
+    Skeleton(model,cell_to_parent_cell,sign_flip)
+  end
+  GridapDistributed.DistributedTriangulation(dtrians,model)
+end
+
+function get_cell_normal_vector(dskeleton::GridapDistributed.DistributedTriangulation)
+  fields=map_parts(dskeleton.trians) do skeleton
+    get_cell_normal_vector(skeleton)
+  end
+  GridapDistributed.DistributedCellField(fields)
+end
+
+function get_cell_owner_normal_vector(dskeleton::GridapDistributed.DistributedTriangulation)
+  fields=map_parts(dskeleton.trians) do skeleton
+    get_cell_owner_normal_vector(skeleton)
+  end
+  GridapDistributed.DistributedCellField(fields)
+end

src/GridapHybrid.jl

@@ -37,4 +37,9 @@ include("HybridLinearSolvers.jl")
 include("GridapAPIExtensions.jl")
 include("GridapTmpModifications.jl")
 
+using GridapDistributed
+using PartitionedArrays
+include("Distributed/Skeleton.jl")
+include("Distributed/HybridAffineFEOperators.jl")
+
 end # module

src/GridapTmpModifications.jl

@@ -302,11 +302,12 @@ function Gridap.Arrays.evaluate!(cache,
 end
 
 function _generate_cell_lface_dofs_from_cell_dofs(glue,
+                                                  cell_grid,
                                                   cell_wise_facets,
                                                   cell_values_field,
                                                   facets_dofs_ids)
   facet_sizes=lazy_map(x->length(x),facets_dofs_ids)
-  cell_wise_facet_sizes=SkeletonVectorFromFacetVector(glue,cell_wise_facets,facet_sizes)
+  cell_wise_facet_sizes=SkeletonVectorFromFacetVector(glue,cell_grid,cell_wise_facets,facet_sizes)
   m=SkeletonVectorFromNonBlockedSkeletonVector(cell_wise_facet_sizes)
   lazy_map(m,cell_values_field,collect(1:length(cell_values_field)))
 end
@@ -335,13 +336,15 @@ function Gridap.FESpaces._change_argument(
         cell_lface_dof_values =
          _generate_cell_lface_dofs_from_cell_dofs(
             trian.glue,
+            trian.grid.parent,
             cell_wise_facets,
             cell_values_field,
             get_cell_dof_ids(Ui))
         Ui_basis = get_fe_basis(Ui)
-        Ui_basis_data = get_data(Ui_basis)
+        Ui_basis_data = Gridap.CellData.get_data(Ui_basis)
         Ui_basis_cell_lface_data = SkeletonVectorFromFacetVector(
           trian.glue,
+          trian.grid.parent,
           cell_wise_facets,
           Ui_basis_data)
         cell_field = lazy_map(linear_combination, cell_lface_dof_values, Ui_basis_cell_lface_data)

src/HybridAffineFEOperators.jl

@@ -1,16 +1,16 @@
 struct HybridAffineFEOperator{TB,TS} <: FEOperator
     weakform::Function
-    trial::MultiFieldFESpace
-    test::MultiFieldFESpace
+    trial
+    test
     bulk_fields::TB
     skeleton_fields::TS
     skeleton_op::AffineFEOperator
 end
 
 function HybridAffineFEOperator(
   weakform::Function,
-  trial::MultiFieldFESpace,
-  test::MultiFieldFESpace,
+  trial,
+  test,
   bulk_fields::TB,
   skeleton_fields::TS) where {TB <: Vector{<:Integer},TS <: Vector{<:Integer}}
 
@@ -27,6 +27,8 @@ function HybridAffineFEOperator(
   # Pair LHS and RHS terms associated to SkeletonTriangulation
     matvec, mat, vec = Gridap.FESpaces._pair_contribution_when_possible(obiform, oliform)
 
+
+
   # Add StaticCondensationMap to matvec terms
     matvec = _add_static_condensation(matvec, bulk_fields, skeleton_fields)
 
@@ -71,7 +73,7 @@ end
 
 function Gridap.FESpaces.solve!(uh, solver::LinearFESolver, op::HybridAffineFEOperator, cache)
   # Solve linear system defined on the skeleton
-    lh = solve(op.skeleton_op)
+    lh = solve(solver,op.skeleton_op)
 
   # Invoke weak form of the hybridizable system
     u = get_trial_fe_basis(op.trial)
@@ -121,25 +123,39 @@ function _compute_hybridizable_from_skeleton_free_dof_values(skeleton_fe_functio
     cell_wise_facets = _get_cell_wise_facets(model)
     cells_around_facets = _get_cells_around_facets(model)
 
+    if (isa(Γ.grid.parent,Gridap.Geometry.GridView))
+      cell_to_parent_cell = Γ.grid.parent.cell_to_parent_cell
+      ifacet_to_facet, facet_to_icell =
+       _find_faces_touched_by_cells(cell_to_parent_cell,
+          cell_wise_facets,
+          cells_around_facets)
+    end
+
     nfields = length(bulk_fields) + length(skeleton_fields)
     m = Gridap.Fields.BlockMap(nfields, skeleton_fields)
     L = Gridap.FESpaces.get_fe_space(skeleton_fe_function)
+    L_cell_dof_ids = get_cell_dof_ids(L)
+    if (isa(Γ.grid.parent,Gridap.Geometry.GridView))
+      L_cell_dof_ids=lazy_map(Reindex(L_cell_dof_ids),ifacet_to_facet)
+    end
+
     lhₑ = lazy_map(m,
                  convert_cell_wise_dofs_array_to_facet_dofs_array(
                  cells_around_facets,
                  cell_wise_facets,
                  lhₖ,
-                 get_cell_dof_ids(L))...)
+                 L_cell_dof_ids)...)
 
     assem = SparseMatrixAssembler(trial_hybridizable, test_hybridizable)
-    lhₑ_dofs = get_cell_dof_ids(trial_hybridizable, get_triangulation(L))
+    Ltrian=get_triangulation(L)
+    if (isa(Γ.grid.parent,Gridap.Geometry.GridView))
+      Ltrian=view(Ltrian,ifacet_to_facet)
+    end
+    lhₑ_dofs = get_cell_dof_ids(trial_hybridizable, Ltrian)
     lhₑ_dofs = lazy_map(m, lhₑ_dofs.args[skeleton_fields]...)
 
-    Ω = trial_hybridizable[first(bulk_fields)]
-    Ω = get_triangulation(Ω)
-
     m = Gridap.Fields.BlockMap(length(bulk_fields), bulk_fields)
-    uhph_dofs = get_cell_dof_ids(trial_hybridizable, Ω)
+    uhph_dofs = get_cell_dof_ids(trial_hybridizable, Γ)
   # This last step is needed as get_cell_dof_ids(...) returns as many blocks
   # as fields in trial, regardless of the FEspaces defined on Ω or not
     uhph_dofs = lazy_map(m, uhph_dofs.args[bulk_fields]...)
@@ -188,40 +204,47 @@ function convert_cell_wise_dofs_array_to_facet_dofs_array(
       cells_around_facets,
       cell_wise_facets,
       cell_dofs::AbstractVector{<:AbstractVector{<:Real}},
-      facet_dofs_ids::AbstractVector{<:AbstractVector{<:Integer}})
-    glue = _generate_glue_among_facet_and_cell_wise_dofs_arrays(
-   cells_around_facets, cell_wise_facets, facet_dofs_ids)
+      facet_dofs_ids::AbstractVector{<:AbstractVector{<:Integer}},
+      ifacet_to_facet::AbstractVector{<:Integer}=IdentityVector(length(facet_dofs_ids)))
+   glue = _generate_glue_among_facet_and_cell_wise_dofs_arrays(
+   cells_around_facets, cell_wise_facets, facet_dofs_ids, ifacet_to_facet)
     k = ExtractFacetDofsFromCellDofsMap(cell_dofs)
     [lazy_map(k, glue)]
 end
 
 function _generate_glue_among_facet_and_cell_wise_dofs_arrays(
   cells_around_facets,
   cell_wise_facets,
-  facet_dof_ids::AbstractVector{<:AbstractVector{<:Integer}})
-
-    c1 = array_cache(cells_around_facets)
-    c2 = array_cache(cell_wise_facets)
-    c3 = array_cache(facet_dof_ids)
-
-    result = Vector{NTuple{3,Int}}(undef, length(facet_dof_ids))
-    current = 1
-    ndofs = 0
-    for facet_gid = 1:length(cells_around_facets)
-        cell_gid = Gridap.Arrays.getindex!(c1, cells_around_facets, facet_gid)[1]
-        current_cell_facets = Gridap.Arrays.getindex!(c2, cell_wise_facets, cell_gid)
-        pos = 1
-        for facet_gid_in_cell in current_cell_facets
-            ndofs = length(Gridap.Arrays.getindex!(c3, facet_dof_ids, facet_gid_in_cell))
-            if (facet_gid == facet_gid_in_cell)
-                break
-            else
-                pos = pos + ndofs
-            end
-        end
-        result[facet_gid] = (cell_gid, pos, ndofs)
-    end
-    result
+  facet_dof_ids::AbstractVector{<:AbstractVector{<:Integer}},
+  ifacet_to_facet::AbstractVector{<:Integer}=IdentityVector(length(facet_dof_ids)))
+
+  @check length(ifacet_to_facet) == length(facet_dof_ids)
+
+  facet_to_ifacet=Dict([facet=>ifacet for (ifacet,facet) in enumerate(ifacet_to_facet)])
+
+  c1 = array_cache(cells_around_facets)
+  c2 = array_cache(cell_wise_facets)
+  c3 = array_cache(facet_dof_ids)
+
+  result = Vector{NTuple{3,Int}}(undef, length(facet_dof_ids))
+  current = 1
+  ndofs = 0
+  for (ifacet,facet_gid) in enumerate(ifacet_to_facet)
+      cell_gid = Gridap.Arrays.getindex!(c1, cells_around_facets, facet_gid)[1]
+      current_cell_facets = Gridap.Arrays.getindex!(c2, cell_wise_facets, cell_gid)
+      pos = 1
+      for facet_gid_in_cell in current_cell_facets
+          ndofs = length(Gridap.Arrays.getindex!(c3, facet_dof_ids, facet_gid_in_cell))
+
+          if (facet_gid == facet_gid_in_cell)
+              break
+          else
+              pos = pos + ndofs
+          end
+      end
+      result[ifacet] = (cell_gid, pos, ndofs)
+  end
+  result
 end
 
 function convert_cell_wise_dofs_array_to_facet_dofs_array(
@@ -362,10 +385,18 @@ function Gridap.FESpaces.get_cell_fe_data(
   tglue::SkeletonGlue) where Dc
     model = tglue.trian.model
     if Dc == num_cell_dims(model)
-        sface_to_data
+      if (isa(tglue.trian.grid.parent,Gridap.Geometry.GridView))
+        cell_to_parent_cell=tglue.trian.grid.parent.cell_to_parent_cell
+        sface_to_data = lazy_map(Reindex(sface_to_data),cell_to_parent_cell)
+      end
+      sface_to_data
     else
         @assert Dc == num_cell_dims(model) - 1
         cell_wise_facets = _get_cell_wise_facets(model)
+        if (isa(tglue.trian.grid.parent,Gridap.Geometry.GridView))
+          cell_to_parent_cell=tglue.trian.grid.parent.cell_to_parent_cell
+          cell_wise_facets = lazy_map(Reindex(cell_wise_facets),cell_to_parent_cell)
+        end
         restrict_facet_dof_ids_to_cell_boundary(cell_wise_facets, sface_to_data)
     end
 end
@@ -441,7 +472,12 @@ end
 function Gridap.FESpaces.get_cell_fe_data(
   fun::typeof(Gridap.FESpaces.get_cell_is_dirichlet),sface_to_data,sglue::FaceToFaceGlue,tglue::SkeletonGlue)
     model = tglue.trian.model
+    grid  = tglue.trian.grid.parent
     cell_wise_facets = _get_cell_wise_facets(model)
+    if (isa(grid,Gridap.Geometry.GridView))
+      cell_to_parent_cell = grid.cell_to_parent_cell
+      cell_wise_facets=lazy_map(Reindex(cell_wise_facets),cell_to_parent_cell)
+    end
     fdofscb = restrict_facet_dof_ids_to_cell_boundary(cell_wise_facets, sface_to_data)
     _generate_cell_is_dirichlet(fdofscb)
 end