Added CartesianDiscreteModels with an arbitrary number of ghost layers

Author: JordiManyer

Project.toml

@@ -5,6 +5,7 @@ version = "0.4.8"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+CircularArrays = "7a955b69-7140-5f4e-a0ed-f168c5e2e749"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
@@ -17,6 +18,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 BlockArrays = "1"
+CircularArrays = "1.4.0"
 FillArrays = "1"
 ForwardDiff = "0.10, 1"
 Gridap = "0.19.1"

src/Geometry.jl

@@ -243,7 +243,11 @@ const DistributedCartesianDiscreteModel{Dc,Dp,A,B,C} =
 function Geometry.CartesianDiscreteModel(
   ranks::AbstractArray{<:Integer}, # Distributed array with the rank IDs
   parts::NTuple{N,<:Integer},      # Number of ranks (parts) in each direction
-  args...;isperiodic=map(i->false,parts),kwargs...) where N
+  args...;
+  ghost = map(i->true,parts),
+  isperiodic = map(i->false,parts),
+  kwargs...
+) where N
 
   desc = CartesianDescriptor(args...;isperiodic=isperiodic,kwargs...)
   nc = desc.partition
@@ -254,10 +258,10 @@ function Geometry.CartesianDiscreteModel(
   @assert N == length(nc) msg
 
   if any(isperiodic)
+    @notimplementedif ghost != map(i->true,parts) 
     _cartesian_model_with_periodic_bcs(ranks,parts,desc)
   else
-    ghost = map(i->true,parts)
-    upartition = uniform_partition(ranks,parts,nc,ghost,isperiodic)
+    upartition = _uniform_partition(ranks,parts,nc,ghost,isperiodic)
     gcids  = CartesianIndices(nc)
     models = map(ranks,upartition) do rank, upartition
       cmin = gcids[first(upartition)]

src/GridapDistributed.jl

@@ -25,6 +25,7 @@ using FillArrays
 using BlockArrays
 using LinearAlgebra
 using ForwardDiff
+using CircularArrays
 
 import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
 import LinearAlgebra: det, tr, cross, dot, ⋅, diag
@@ -43,6 +44,8 @@ export with_ghost, no_ghost
 
 export redistribute
 
+include("PArraysExtras.jl")
+
 include("BlockPartitionedArrays.jl")
 
 include("Algebra.jl")

src/PArraysExtras.jl

@@ -0,0 +1,84 @@
+
+# This is copied from PArrays v0.5, which is not yet supported. It is necessary to get 
+# more than one layers of ghosts.
+
+function _uniform_partition(rank,np,n,args...)
+  @assert prod(np) == length(rank)
+  indices = map(rank) do rank
+    _block_with_constant_size(rank,np,n,args...)
+  end
+  if length(args) == 0
+    map(indices) do indices
+      cache = assembly_cache(indices)
+      copy!(cache,empty_assembly_cache())
+    end
+  else
+    assembly_neighbors(indices;symmetric=true)
+  end
+  indices
+end
+
+function _block_with_constant_size(rank,np,n,ghost,periodic=map(i->false,ghost))
+  N = length(n)
+  p = CartesianIndices(np)[rank]
+  own_ranges = map(_local_range,Tuple(p),np,n)
+  local_ranges = map(_local_range,Tuple(p),np,n,ghost,periodic)
+  owners = map(Tuple(p), np, n, local_ranges) do p, np, n, lr
+    myowners = zeros(Int32,length(lr))
+    i = 1
+    for p in Iterators.cycle(1:np)
+      plr = _local_range(p, np, n)
+      while mod(lr[i]-1, n)+1 in plr
+        myowners[i] = p
+        (i += 1) > length(myowners) && return myowners
+      end
+    end
+  end
+  n_local = prod(map(length, local_ranges))
+  n_own = prod(map(length, own_ranges))
+  n_ghost = n_local - n_own
+
+  ghost_to_global = zeros(Int,n_ghost)
+  ghost_to_owner = zeros(Int32,n_ghost)
+  perm = zeros(Int32,n_local)
+  i_ghost = 0
+  i_own = 0
+
+  cis = CartesianIndices(map(length,local_ranges))
+  lis = CircularArray(LinearIndices(n))
+  local_cis = CartesianIndices(local_ranges)
+  owner_lis = LinearIndices(np)
+  for (i,ci) in enumerate(cis)
+    flags = map(Tuple(ci), own_ranges, local_ranges) do i, or, lr
+      i in (or .- first(lr) .+ 1)
+    end
+    if !all(flags)
+      i_ghost += 1
+      ghost_to_global[i_ghost] = lis[local_cis[i]]
+      o = map(getindex,owners,Tuple(ci))
+      o_ci = CartesianIndex(o)
+      ghost_to_owner[i_ghost] = owner_lis[o_ci]
+      perm[i] = i_ghost + n_own
+    else
+      i_own += 1
+      perm[i] = i_own
+    end
+  end
+  ghostids = GhostIndices(prod(n),ghost_to_global,ghost_to_owner)
+  ids = PartitionedArrays.LocalIndicesWithConstantBlockSize(p,np,n,ghostids)
+  PartitionedArrays.PermutedLocalIndices(ids,perm)
+end
+
+function _local_range(p,np,n,ghost=false,periodic=false)
+  l, rem = divrem(n, np)
+  offset = l * (p-1)
+  if rem >= (np-p+1)
+    l += 1
+    offset += p - (np-rem) - 1
+  end
+  start = 1+offset-ghost
+  stop = l+offset+ghost
+
+  periodic && return start:stop
+  return max(1, start):min(n,stop)
+end

test/ExtraGhostLayers.jl

@@ -0,0 +1,36 @@
+
+using Gridap
+using GridapDistributed
+using PartitionedArrays
+
+
+parts = (2,2)
+ranks = with_debug() do distribute
+  distribute(LinearIndices((prod(parts),)))
+end
+
+nc = (4,4)
+model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),nc;ghost=(2,2))
+
+num_cells(model)
+map(num_cells,local_views(model))
+
+reffe = ReferenceFE(lagrangian,Float64,1)
+V = FESpace(model,reffe)
+
+cell_dof_ids = map(get_cell_dof_ids,local_views(V))
+
+gids = get_free_dof_ids(V)
+own_to_local(gids)
+
+Ω = Triangulation(model)
+map(num_cells,local_views(Ω))
+
+dΩ = Measure(Ω,2)
+
+f = 1.0
+a(u,v) = ∫(∇(u)⋅∇(v))dΩ
+l(v) = ∫(f⋅v)dΩ
+
+
+