gridap / gridap/GridapDistributed.jl

Potential issue when interpolating DistributedCelllFields that are only defined on own cells

Open
#184 1 comment 0 reactions 2 assignees Claimed by @JordiManyer View on GitHub
Dominant language
Julia
Stars
148
Forks
28
PR merge metrics
No merged PRs in 30d

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).

```julia
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
```

![Image](https://github.com/user-attachments/assets/4b4301dc-6cd8-4d7f-8a71-de58842ab3db)

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.