DistributedStokesTests.jl and MPIPETScDistributedStokesTests.jl working

Author: amartinhuertas

src/DistributedAssemblers.jl

@@ -496,17 +496,17 @@ end
 function Gridap.Geometry.Triangulation(strategy::Type{<:AssemblyStrategy},
   model::DistributedDiscreteModel,args_triangulation...)#;kwargs_triangulation...)
   DistributedData(model) do part, (model,gids)
-    Triangulation(strategy,model,part,gids,args_triangulation...)#;kwargs_triangulation)
+    Triangulation(strategy,part,gids,model,args_triangulation...)#;kwargs_triangulation)
   end
 end
 
 function Gridap.Geometry.Triangulation(strategy::AssemblyStrategy,
   part,gids,model::DiscreteModel,args_triangulation...)#;kwargs_triangulation...)
-  Triangulation(typeof(strategy),model,part,gids,args_triangulation...)#;kwargs_triangulation)
+  Triangulation(typeof(strategy),part,gids,model,args_triangulation...)#;kwargs_triangulation)
 end
 
 function Gridap.Geometry.Triangulation(strategy::Type{<:AssemblyStrategy},
-  model::DiscreteModel,part,gids,args_triangulation...)#;kwargs_triangulation...)
+  part,gids,model::DiscreteModel,args_triangulation...)#;kwargs_triangulation...)
   portion=proc_local_triangulation_portion_type(strategy)
   Triangulation(portion,part,gids,model,args_triangulation...)#;kwargs_triangulation)
 end
@@ -518,8 +518,32 @@ function Gridap.Geometry.Triangulation(strategy::DistributedAssemblyStrategy,
   end
 end
 
-# function Gridap.Geometry.Triangulation(strategy::AssemblyStrategy,
-#   model::DiscreteModel,args_triangulation...)#;kwargs_triangulation...)
-#   portion_type=get_proc_local_triangulation_portion_type(typeof(strategy))
-#   Triangulation(portion_type, model, args_triangulation...)#;kwargs_triangulation...)
-# end
+function Gridap.Geometry.BoundaryTriangulation(model::DistributedDiscreteModel;kwargs_triangulation...)
+  AS=default_assembly_strategy_type(get_comm(model))
+  BoundaryTriangulation(AS, model; kwargs_triangulation...)
+end
+
+function Gridap.Geometry.BoundaryTriangulation(strategy::Type{<:AssemblyStrategy},
+  model::DistributedDiscreteModel;kwargs_triangulation...)
+  DistributedData(model) do part, (model,gids)
+    BoundaryTriangulation(strategy,part,gids,model;kwargs_triangulation...)
+  end
+end
+
+function Gridap.Geometry.BoundaryTriangulation(strategy::AssemblyStrategy,
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  BoundaryTriangulation(typeof(strategy),part,gids,model;kwargs_triangulation...)
+end
+
+function Gridap.Geometry.BoundaryTriangulation(strategy::Type{<:AssemblyStrategy},
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  portion=proc_local_triangulation_portion_type(strategy)
+  BoundaryTriangulation(portion,part,gids,model;kwargs_triangulation...)
+end
+
+function Gridap.Geometry.BoundaryTriangulation(strategy::DistributedAssemblyStrategy,
+  model::DistributedDiscreteModel;kwargs_triangulation...)
+  DistributedData(model,strategy) do part, (model,gids), strategy
+     BoundaryTriangulation(strategy,part,gids,model;kwargs_triangulation...)
+  end
+end

src/DistributedDiscreteModels.jl

@@ -24,4 +24,33 @@ function Gridap.writevtk(model::DistributedDiscreteModel,filebase::String)
 
 end
 
+function Gridap.Geometry.get_face_labeling(model::DistributedDiscreteModel)
+  DistributedData(model) do part, (model,gids)
+     get_face_labeling(model)
+  end
+end
+
+"""
+    add_tag_from_tags!(lab::FaceLabeling, name::String, tags::Vector{Int})
+    add_tag_from_tags!(lab::FaceLabeling, name::String, tags::Vector{String})
+    add_tag_from_tags!(lab::FaceLabeling, name::String, tag::Int)
+    add_tag_from_tags!(lab::FaceLabeling, name::String, tag::String)
+"""
+function Gridap.Geometry.add_tag_from_tags!(lab::DistributedData{<:FaceLabeling}, name::String, tags::Vector{Int})
+  do_on_parts(lab) do part, lab
+    add_tag_from_tags!(lab, name, tags)
+  end
+end
 
+function Gridap.Geometry.add_tag_from_tags!(
+  labels::DistributedData{<:FaceLabeling}, name::String, names::Vector{String})
+  @notimplemented
+end
+
+function Gridap.Geometry.add_tag_from_tags!(labels::DistributedData{<:FaceLabeling}, name::String, tag::Int)
+  @notimplemented
+end
+
+function Gridap.Geometry.add_tag_from_tags!(labels::DistributedData{<:FaceLabeling}, name::String, tag::String)
+  @notimplemented
+end

src/DistributedFESpaces.jl

@@ -388,12 +388,8 @@ struct DistributedFEFunction <: FEFunction
     space::DistributedFESpace
 end
 
-# TO-DO Gridap.FESpaces.FEFunctionStyle(::Type{DistributedFEFunction}) = Val{true}()
-
 get_distributed_data(u::DistributedFEFunction) = u.funs
 
 Gridap.FESpaces.get_free_dof_values(a::DistributedFEFunction) = a.vals
 
 Gridap.FESpaces.get_fe_space(a::DistributedFEFunction) = a.space
-
-# TO-DO Gridap.FESpaces.is_a_fe_function(a::DistributedFEFunction) = true

src/DistributedTriangulations.jl

@@ -53,20 +53,34 @@ end
 
 #TO-DO: what does it mean OwnedCells and OwnedAndGhostCells for a BoundaryTriangulation?
 #       Perhaps we should use a different type, with a more intention revealing name!
+function Gridap.Geometry.BoundaryTriangulation(portion::Type{<:ProcLocalTriangulationPortion},
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  @abstractmethod
+end
+
+function Gridap.Geometry.BoundaryTriangulation(portion::Type{<:ProcLocalTriangulationPortion},
+  model::DistributedDiscreteModel;kwargs_triangulation...)
+  DistributedData(model) do part, (model,gids)
+    BoundaryTriangulation(portion,part,gids,model;kwargs_triangulation...)
+  end
+end
+
 function Gridap.Geometry.BoundaryTriangulation(portion::Type{OwnedCells},
-  model::DiscreteModel,args_triangulation...)#;kwargs_triangulation...)
-  @notimplemented
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  trian=BoundaryTriangulation(model;kwargs_triangulation...)
+  filter_cells_when_needed(portion,trian,part,gids)
 end
 
 function Gridap.Geometry.BoundaryTriangulation(portion::Type{OwnedAndGhostCells},
-  model::DistributedDiscreteModel,args_triangulation...)#;kwargs_triangulation...)
-  @notimplemented
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  trian=BoundaryTriangulation(model;kwargs_triangulation...)
+  filter_cells_when_needed(portion,trian)
 end
 
 #TO-DO: what does it mean OwnedCells and OwnedAndGhostCells for a SkeletonTriangulation?
 #       Perhaps we should use a different type, with a more intention revealing name!
 function Gridap.Geometry.SkeletonTriangulation(portion::Type{OwnedCells},
-  model::DistributedDiscreteModel,args_triangulation...;kwargs_triangulation...)
+  model::DistributedDiscreteModel;kwargs_triangulation...)
   @notimplemented
 end
 

src/GridapHighLevelAPI.jl

@@ -117,6 +117,28 @@ for T in (DistributedData{<:CellField},DistributedFEFunction)
   end
  end
 
+function (+)(a::DistributedData{<:Gridap.CellData.DomainContribution},
+             b::DistributedData{<:Gridap.CellData.DomainContribution})
+   DistributedData(a,b) do part, a, b
+     a+b
+   end
+end
+
+function (-)(a::DistributedData{<:Gridap.CellData.DomainContribution},
+             b::DistributedData{<:Gridap.CellData.DomainContribution})
+   DistributedData(a,b) do part, a, b
+     a-b
+   end
+end
+
+function (*)(a::Number,b::DistributedData{<:Gridap.CellData.DomainContribution})
+  DistributedData(b) do part, b
+    a*b
+  end
+end
+
+(*)(a::DistributedData{<:Gridap.CellData.DomainContribution},b::Number) = b*a
+
  function (a::typeof(gradient))(x::DistributedData{<:CellField})
    DistributedData(x) do part, x
      a(x)
@@ -128,3 +150,40 @@ for T in (DistributedData{<:CellField},DistributedFEFunction)
     a(x)
   end
 end
+
+function dot(::typeof(∇),f::DistributedData{<:CellField})
+  DistributedData(f) do part, f
+    divergence(f)
+  end
+end
+
+
+function Base.iterate(a::DistributedData{<:MultiFieldCellField})
+  if _num_fields(a)==0
+    return nothing
+  end
+  sf=_get_field(a,1)
+  (sf,2)
+end
+
+function Base.iterate(a::DistributedData{<:MultiFieldCellField},state)
+  if state > _num_fields(a)
+     return nothing
+  end
+  sf=_get_field(a,state)
+  (sf,state+1)
+end
+
+function _num_fields(a::DistributedData{<:MultiFieldCellField})
+  num_fields=0
+  do_on_parts(a) do part, a
+    num_fields=length(a.single_fields)
+  end
+  num_fields
+end
+
+function _get_field(a::DistributedData{<:MultiFieldCellField},field_id)
+  DistributedData(a) do part, a
+     a.single_fields[field_id]
+  end
+end

src/MultiFieldDistributedFESpaces.jl

@@ -207,8 +207,6 @@ struct MultiFieldDistributedFEFunction <: FEFunction
     space::MultiFieldDistributedFESpace
 end
 
-# TO-DO Gridap.FESpaces.FEFunctionStyle(::Type{MultiFieldDistributedFEFunction}) = Val{true}()
-
 get_distributed_data(u::MultiFieldDistributedFEFunction) =
      get_distributed_data(u.multifield_fe_function)
 
@@ -217,4 +215,6 @@ Gridap.FESpaces.get_free_dof_values(a::MultiFieldDistributedFEFunction) =
 
 Gridap.FESpaces.get_fe_space(a::MultiFieldDistributedFEFunction) = a.space
 
-# TO-DO Gridap.FESpaces.is_a_fe_function(a::MultiFieldDistributedFEFunction) = true
+Base.iterate(m::MultiFieldDistributedFEFunction)=iterate(m.single_fe_functions)
+Base.iterate(m::MultiFieldDistributedFEFunction,state) = iterate(m.single_fe_functions,state)
+Base.getindex(m::MultiFieldDistributedFEFunction,field_id::Integer) = m.single_fe_functions[field_id]

test/DistributedStokesTests.jl

@@ -15,13 +15,6 @@ function Gridap.FESpaces.num_dirichlet_dofs(f::Gridap.MultiField.MultiFieldFESpa
 end
 
 function run(comm,subdomains)
-  # 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
   ux(x)=2*x[1]*x[2]
   uy(x)=-x[2]^2
@@ -37,29 +30,23 @@ function run(comm,subdomains)
   cells = (4, 4)
   model = CartesianDiscreteModel(comm, subdomains, domain, cells)
 
-  # Define Dirichlet and Neumann boundaries for local models
-  do_on_parts(model) do part, (model,gids)
-    labels = get_face_labeling(model)
-    add_tag_from_tags!(labels,"diri0",[1,2,3,4,6,7,8])
-    add_tag_from_tags!(labels,"diri1",[1,2,3,4,6,7,8])
-    add_tag_from_tags!(labels,"neumann",[5])
-  end
-
+  labels = get_face_labeling(model)
+  add_tag_from_tags!(labels,"diri0",[1,2,3,4,6,7,8])
+  add_tag_from_tags!(labels,"diri1",[1,2,3,4,6,7,8])
+  add_tag_from_tags!(labels,"neumann",[5])
 
   # FE Spaces
   order = 2
   reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
   reffeₚ = Gridap.ReferenceFEs.LagrangianRefFE(Float64,QUAD,order-1;space=:P)
   V = FESpace(
-        vector_type,
         model=model,
         reffe=reffeᵤ,
         conformity=:H1,
         dirichlet_tags=["diri0","diri1"],
       )
 
   Q = FESpace(
-       vector_type,
        model=model,
        reffe=reffeₚ,
        conformity=:L2)
@@ -71,53 +58,30 @@ function run(comm,subdomains)
   P=TrialFESpace(Q)
   X=MultiFieldFESpace(Y,[U,P])
 
-  strategy = OwnedAndGhostCellsAssemblyStrategy(V,MapDoFsTypeGlobal())
+  trian=Triangulation(model)
+  degree = 2*(order+1)
+  dΩ = Measure(trian,degree)
 
-  function a(x,y)
-    DistributedData(x,y,ddΩ,ddΓ) do part, xl, yl, dΩ, dΓ
-        ul,pl=xl
-        vl,ql=yl
-        ∫( ∇(vl)⊙∇(ul) - (∇⋅vl)*pl + ql*(∇⋅ul) )dΩ
-    end
-  end
+  btrian=BoundaryTriangulation(model;tags="neumann")
+  dΓ = Measure(btrian,degree)
 
-  function l(y)
-    DistributedData(y,ddΩ,ddΓ) do part, yl, dΩ, dΓ
-      vl,_=yl
-      ∫( vl⋅f )dΩ + ∫( vl⋅s )dΓ
-    end
+  function a((u,p),(v,q))
+    ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )dΩ
   end
 
-   # Assembler
-   assem = SparseMatrixAssembler(matrix_type, vector_type, X, Y, strategy)
-
-   function setup_dΩ(part,(model,gids),strategy)
-     trian = Triangulation(strategy,model)
-     degree = 2*(order+1)
-     Measure(trian,degree)
-   end
-   ddΩ = DistributedData(setup_dΩ,model,strategy)
-
-   function setup_dΓ(part,(model,gids),strategy)
-     btrian=BoundaryTriangulation(strategy,model;tags="neumann")
-     degree = 2*(order+1)
-     Measure(btrian,degree)
-   end
-   ddΓ = DistributedData(setup_dΓ,model,strategy)
+  function l((v,q))
+    ∫( v⋅f )dΩ + ∫( v⋅s )dΓ
+  end
 
   # #  # FE solution
-  op = AffineFEOperator(a,l,X,Y,assem)
+  op = AffineFEOperator(a,l,X,Y)
   xh = solve(op)
+  uh,_ = xh
 
-  sums = DistributedData(model, xh) do part, (model, gids), xh
-    trian = Triangulation(model)
-    owned_trian = remove_ghost_cells(trian, part, gids)
-    dΩ = Measure(owned_trian, 2*order)
-    uh,_ = xh
-    e = u-uh
-    sum(∫(e⋅e)dΩ)
-  end
-  e_l2 = sum(gather(sums))
+  trian=Triangulation(OwnedCells,model)
+  dΩ=Measure(trian,2*order)
+  e = u-uh
+  e_l2 = sum(∫(e⋅e)dΩ)
   tol = 1.0e-9
   println("$(e_l2) < $(tol)")
   @test e_l2 < tol
@@ -128,6 +92,6 @@ function f(comm)
   run(comm,subdomains)
 end
 
-#SequentialCommunicator(f,subdomains)
+SequentialCommunicator(f,subdomains)
 
 end # module