DistributedPoissonDGTests.jl and MPIPETScDistributedPoissonDGTests.jl

working with high level API

Author: amartinhuertas

src/DistributedAssemblers.jl

@@ -544,6 +544,36 @@ 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...)
+    BoundaryTriangulation(strategy,part,gids,model;kwargs_triangulation...)
+  end
+end
+
+function Gridap.Geometry.SkeletonTriangulation(model::DistributedDiscreteModel;kwargs_triangulation...)
+  AS=default_assembly_strategy_type(get_comm(model))
+  SkeletonTriangulation(AS, model; kwargs_triangulation...)
+end
+
+function Gridap.Geometry.SkeletonTriangulation(strategy::Type{<:AssemblyStrategy},
+  model::DistributedDiscreteModel;kwargs_triangulation...)
+  DistributedData(model) do part, (model,gids)
+    SkeletonTriangulation(strategy,part,gids,model;kwargs_triangulation...)
+  end
+end
+
+function Gridap.Geometry.SkeletonTriangulation(strategy::AssemblyStrategy,
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  SkeletonTriangulation(typeof(strategy),part,gids,model;kwargs_triangulation...)
+end
+
+function Gridap.Geometry.SkeletonTriangulation(strategy::Type{<:AssemblyStrategy},
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  portion=proc_local_triangulation_portion_type(strategy)
+  SkeletonTriangulation(portion,part,gids,model;kwargs_triangulation...)
+end
+
+function Gridap.Geometry.SkeletonTriangulation(strategy::DistributedAssemblyStrategy,
+  model::DistributedDiscreteModel;kwargs_triangulation...)
+  DistributedData(model,strategy) do part, (model,gids), strategy
+     SkeletonTriangulation(strategy,part,gids,model;kwargs_triangulation...)
   end
 end

src/DistributedTriangulations.jl

@@ -79,14 +79,28 @@ 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},
+function Gridap.Geometry.SkeletonTriangulation(portion::Type{<:ProcLocalTriangulationPortion},
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  @abstractmethod
+end
+
+function Gridap.Geometry.SkeletonTriangulation(portion::Type{<:ProcLocalTriangulationPortion},
   model::DistributedDiscreteModel;kwargs_triangulation...)
-  @notimplemented
+  DistributedData(model) do part, (model,gids)
+    SkeletonTriangulation(portion,part,gids,model;kwargs_triangulation...)
+  end
+end
+
+function Gridap.Geometry.SkeletonTriangulation(portion::Type{OwnedCells},
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  trian=SkeletonTriangulation(model;kwargs_triangulation...)
+  filter_cells_when_needed(portion,trian,part,gids)
 end
 
 function Gridap.Geometry.SkeletonTriangulation(portion::Type{OwnedAndGhostCells},
-  model::DistributedDiscreteModel,args_triangulation...;kwargs_triangulation...)
-  @notimplemented
+  part,gids,model::DiscreteModel;kwargs_triangulation...)
+  trian=SkeletonTriangulation(model;kwargs_triangulation...)
+  filter_cells_when_needed(portion,trian)
 end
 
 function remove_ghost_cells(trian::Triangulation, part::Integer, gids::IndexSet)
@@ -125,3 +139,9 @@ end
 function include_ghost_cells(trian::RestrictedTriangulation)
     trian.oldtrian
 end
+
+function Gridap.CellData.get_normal_vector(trian::DistributedData{<:Triangulation})
+  DistributedData(trian) do part, trian
+    Gridap.CellData.get_normal_vector(trian)
+  end
+end

src/GridapHighLevelAPI.jl

@@ -187,3 +187,36 @@ function _get_field(a::DistributedData{<:MultiFieldCellField},field_id)
      a.single_fields[field_id]
   end
 end
+
+function Gridap.CellData.jump(a::DistributedData{<:CellField})
+  DistributedData(a) do part, a
+     jump(a)
+  end
+end
+
+function Gridap.CellData.jump(a::DistributedData{<:SkeletonPair{<:CellField}})
+  DistributedData(a) do part, a
+     jump(a)
+  end
+end
+
+function Gridap.CellData.mean(a::DistributedData{<:CellField})
+  DistributedData(a) do part, a
+     mean(a)
+  end
+end
+
+for op in (:outer,:*,:dot)
+  @eval begin
+    function ($op)(a::DistributedData{<:CellField},b::DistributedData{<:SkeletonPair{<:CellField}})
+      DistributedData(a,b) do part, a, b
+        Operation($op)(a,b)
+      end
+    end
+    function ($op)(a::DistributedData{<:SkeletonPair{<:CellField}},b::DistributedData{<:CellField})
+      DistributedData(a,b) do part, a, b
+        Operation($op)(a,b)
+      end
+    end
+  end
+end

test/DistributedPoissonDGTests.jl

@@ -47,41 +47,26 @@ function setup_fe_spaces(model)
 end
 
 function run(comm,model,U,V,strategy)
+  trian=Triangulation(strategy,model)
+  dΩ=Measure(trian,degree)
 
-  function setup_dΩ(part,(model,gids),strategy)
-    trian = Triangulation(strategy,model)
-    Measure(trian,degree)
-  end
-  ddΩ = DistributedData(setup_dΩ,model,strategy)
+  btrian=BoundaryTriangulation(strategy,model)
+  dΓ = Measure(btrian,degree)
 
-  function setup_dΓ(part,(model,gids),strategy)
-    trian = BoundaryTriangulation(strategy,model)
-    Measure(trian,degree)
-  end
-  ddΓ = DistributedData(setup_dΓ,model,strategy)
-
-  function setup_dΛ(part,(model,gids),strategy)
-    trian = SkeletonTriangulation(strategy,model)
-    Measure(trian,degree)
-  end
-  ddΛ = DistributedData(setup_dΛ,model,strategy)
+  strian=SkeletonTriangulation(strategy,model)
+  dΛ = Measure(strian,degree)
 
+  n_Γ = get_normal_vector(btrian)
+  n_Λ = get_normal_vector(strian)
 
   function a(u,v)
-    DistributedData(u,v,ddΩ,ddΓ,ddΛ) do part, ul, vl, dΩ, dΓ, dΛ
-      n_Γ = get_normal_vector(dΓ.quad.trian)
-      n_Λ = get_normal_vector(dΛ.quad.trian)
-      ∫( ∇(vl)⋅∇(ul) )*dΩ +
-      ∫( (γ/h)*vl*ul - vl*(n_Γ⋅∇(ul)) - (n_Γ⋅∇(vl))*ul )*dΓ +
-      ∫( (γ/h)*jump(vl*n_Λ)⋅jump(ul*n_Λ) - jump(vl*n_Λ)⋅mean(∇(ul)) -  mean(∇(vl))⋅jump(ul*n_Λ) )*dΛ
-    end
+      ∫( ∇(v)⋅∇(u) )*dΩ +
+      ∫( (γ/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )*dΓ +
+      ∫( (γ/h)*jump(v*n_Λ)⋅jump(u*n_Λ) - jump(v*n_Λ)⋅mean(∇(u)) -  mean(∇(v))⋅jump(u*n_Λ) )*dΛ
   end
   function l(v)
-    DistributedData(v,ddΩ,ddΓ) do part, vl, dΩ, dΓ
-      n_Γ = get_normal_vector(dΓ.quad.trian)
-      ∫( vl*f )*dΩ +
-      ∫( (γ/h)*vl*ud - (n_Γ⋅∇(vl))*ud )*dΓ
-    end
+      ∫( v*f )*dΩ +
+      ∫( (γ/h)*v*ud - (n_Γ⋅∇(v))*ud )*dΓ
   end
 
   # Assembly
@@ -91,21 +76,13 @@ function run(comm,model,U,V,strategy)
   # 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)
-    dΩ = Measure(owned_trian, degree)
-    e = u - uh
-    l2(u) = ∫( u⊙u )*dΩ
-    h1(u) = ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ
-    e_l2 = sum(l2(e))
-    e_h1 = sum(h1(e))
-    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))
+  trian = Triangulation(OwnedCells,model)
+  dΩ = Measure(trian,degree)
+  e = u - uh
+  l2(u) = ∫( u⊙u )*dΩ
+  h1(u) = ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ
+  e_l2 = sum(l2(e))
+  e_h1 = sum(h1(e))
   tol = 1.0e-9
   println("$(e_l2) < $(tol)")
   println("$(e_h1) < $(tol)")

test/MPIPETScTests/MPIPETScDistributedPoissonDGTests.jl

@@ -48,41 +48,26 @@ function setup_fe_spaces(model)
 end
 
 function run(comm,model,U,V,strategy)
+  trian=Triangulation(strategy,model)
+  dΩ=Measure(trian,degree)
 
-  function setup_dΩ(part,(model,gids),strategy)
-    trian = Triangulation(strategy,model)
-    Measure(trian,degree)
-  end
-  ddΩ = DistributedData(setup_dΩ,model,strategy)
+  btrian=BoundaryTriangulation(strategy,model)
+  dΓ = Measure(btrian,degree)
 
-  function setup_dΓ(part,(model,gids),strategy)
-    trian = BoundaryTriangulation(strategy,model)
-    Measure(trian,degree)
-  end
-  ddΓ = DistributedData(setup_dΓ,model,strategy)
-
-  function setup_dΛ(part,(model,gids),strategy)
-    trian = SkeletonTriangulation(strategy,model)
-    Measure(trian,degree)
-  end
-  ddΛ = DistributedData(setup_dΛ,model,strategy)
+  strian=SkeletonTriangulation(strategy,model)
+  dΛ = Measure(strian,degree)
 
+  n_Γ = get_normal_vector(btrian)
+  n_Λ = get_normal_vector(strian)
 
   function a(u,v)
-    DistributedData(u,v,ddΩ,ddΓ,ddΛ) do part, ul, vl, dΩ, dΓ, dΛ
-      n_Γ = get_normal_vector(dΓ.quad.trian)
-      n_Λ = get_normal_vector(dΛ.quad.trian)
-      ∫( ∇(vl)⋅∇(ul) )*dΩ +
-      ∫( (γ/h)*vl*ul - vl*(n_Γ⋅∇(ul)) - (n_Γ⋅∇(vl))*ul )*dΓ +
-      ∫( (γ/h)*jump(vl*n_Λ)⋅jump(ul*n_Λ) - jump(vl*n_Λ)⋅mean(∇(ul)) -  mean(∇(vl))⋅jump(ul*n_Λ) )*dΛ
-    end
+      ∫( ∇(v)⋅∇(u) )*dΩ +
+      ∫( (γ/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )*dΓ +
+      ∫( (γ/h)*jump(v*n_Λ)⋅jump(u*n_Λ) - jump(v*n_Λ)⋅mean(∇(u)) -  mean(∇(v))⋅jump(u*n_Λ) )*dΛ
   end
   function l(v)
-    DistributedData(v,ddΩ,ddΓ) do part, vl, dΩ, dΓ
-      n_Γ = get_normal_vector(dΓ.quad.trian)
-      ∫( vl*f )*dΩ +
-      ∫( (γ/h)*vl*ud - (n_Γ⋅∇(vl))*ud )*dΓ
-    end
+      ∫( v*f )*dΩ +
+      ∫( (γ/h)*v*ud - (n_Γ⋅∇(v))*ud )*dΓ
   end
 
   # Assembly
@@ -101,21 +86,13 @@ function run(comm,model,U,V,strategy)
   fels = LinearFESolver(ls)
   uh = solve(fels, 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)
-    dΩ = Measure(owned_trian, degree)
-    e = u - uh
-    l2(u) = ∫( u⊙u )*dΩ
-    h1(u) = ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ
-    e_l2 = sum(l2(e))
-    e_h1 = sum(h1(e))
-    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))
+  trian = Triangulation(OwnedCells,model)
+  dΩ = Measure(trian,degree)
+  e = u - uh
+  l2(u) = ∫( u⊙u )*dΩ
+  h1(u) = ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ
+  e_l2 = sum(l2(e))
+  e_h1 = sum(h1(e))
   tol = 1.0e-9
   @test e_l2 < tol
   @test e_h1 < tol