Potential issue when interpolating DistributedCelllFields that are only defined on own cells #184
Description
Documenting here an issue that has arisen in GridapP4est.jl/GridapGeosciences.jl when dealing with periodic distributed meshes such that owned cells are locally numbered before ghost cells.
The mwe below illustrates (although it does not reproduce) the issue (as it does not actually use a periodic distributed mesh with the features described in the previous paragraph). This mwe is to be executed with 2 MPI tasks. The triangulation Ω only includes owned cells. As a consequence, the DistributedCellField dcf is only defined on owned cells. When we interpolate the local portions of such DistributedCellField on the FESpace, the underlying code under the interpolate call (extracted here in the mwe for illustration purposes), extends the cell-wise array resulting from the cell field interpolation with zeros on the ghost cells. This particularly happens in the CellField(cf,trian,DomainStyle(dof_basis)) call. Here, the triangulation associated to the FESpace has both owned and ghost cells, while the cell field only has owned cells.
Extending with zeros the ghost cells is dangerous, as the algorithm that gathers the cell-wise DoFs into the DoF-wise array, performs a sweep starting from the first to the last cell. If the mesh is periodic, you may have (1) ghost DoFs on ghost cells that are also on owned cells, and (2) this may happen from the perspective of all processors that share such DoFs. As a consequence the consistent! call in GridapDistributed that typically fixes that with non-periodic meshes, cannot fix that with periodic meshes.
The picture below the mwe illustrates the particular scenario in which the issue is arising. The periodic mesh is distributed among 2 processors. The picture shows the local numbering of DoFs on the first processor (rank 1). The ghost DoFs on owned cells in rank 1 (marked with a cross x in the picture), end up having the value zero after consistent! (because they were overridden with zeros on both processors).
The issue can by by-passed by creating a DistributedCellField defined on a triangulation that has both owned and ghost cells. But is there a better way to handle this such that the user does not have to bother about these details? (I am thinking about this).
using Gridap
using GridapDistributed
using PartitionedArrays
using MPI
using FillArrays
MPI.Init()
ranks = distribute_with_mpi(LinearIndices((prod(MPI.Comm_size(MPI.COMM_WORLD)),)))
model = CartesianDiscreteModel(ranks, (2,1), (0,1,0,1), (4,4))
reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(model,reffe,conformity=:H1)
Ω = Triangulation(model)
fields = map(local_views(Ω)) do Ω
println(num_cells(Ω))
Gridap.CellData.GenericCellField(Fill(Gridap.Fields.ConstantField(1.0),num_cells(Ω)), Ω, ReferenceDomain())
end
dcf = GridapDistributed.DistributedCellField(fields, Ω)
Udofs = get_free_dof_values(interpolate(dcf, V))
map(ranks, local_views(V), local_views(dcf)) do rank, V, cf
dof_basis=Gridap.FESpaces.get_fe_dof_basis(V)
trian = get_triangulation(V)
cfp = CellField(cf,trian,DomainStyle(dof_basis))
x=dof_basis(cfp)
if rank==1
println("Rank $(MPI.Comm_rank(MPI.COMM_WORLD)) x at dofs: ", x); print("\n");
end
end 