Generating PVD files when the Triangulation is not present on all Ranks

[shreyas02]

With the current implementation, we can not generate the paraview files when the triangulations are not present on all MPI ranks. In issue #70, this has been written as a #TODO task for taking nothing.

Since I needed this functionality for my use case, I tried to fix it myself, but I had to use functions from Mpi.jl as well. Below are the changes I made to make it work.

struct DistributedPvd{T<:AbstractArray}
  pvds::T
  main_rank::Int
end

function Visualization.createpvd(trian::DistributedTriangulation,parts::AbstractArray,args...;kwargs...)
  main_rank = -1
  map(local_views(trian), parts) do LocTrian, part
    flag = part * !(num_cells(LocTrian) == 0)
    global_flag = Vector{Int}(undef, MPI.Comm_size(MPI.COMM_WORLD))
    MPI.Allgather!([flag], global_flag, MPI.COMM_WORLD)
    main_rank = maximum(global_flag)
  end
  pvds = map_main(parts,main = main_rank) do part
    p = paraview_collection(args...;kwargs...)
  end
  DistributedPvd(pvds,main_rank)
end

function Visualization.createpvd(f,trian::DistributedTriangulation,parts::AbstractArray,args...;kwargs...)
  pvd = createpvd(trian,parts,args...;kwargs...)
  try
    f(pvd)
  finally
    savepvd(pvd)
  end
end

function Visualization.createpvd(parts::AbstractArray,args...;kwargs...)
  pvds = map_main(parts, main = 1) do part
    p = paraview_collection(args...;kwargs...)
  end
  DistributedPvd(pvds,1)
end

function Visualization.createpvd(f,parts::AbstractArray,args...;kwargs...)
  pvd = createpvd(parts,args...;kwargs...)
  try
    f(pvd)
  finally
    savepvd(pvd)
  end
end

function Visualization.savepvd(pvd::DistributedPvd)
  map_main(pvd.pvds, main = pvd.main_rank) do pvd
    vtk_save(pvd)
  end
end

function Base.setindex!(pvd::DistributedPvd,pvtk::AbstractArray,time::Real)

  map(pvtk) do pvtk
    if !isnothing(pvtk)
      vtk_save(pvtk)
    end  
  end

  map_main(pvtk,pvd.pvds,main = pvd.main_rank) do pvtk,pvd
    pvd[time] = pvtk
  end
end

[JordiManyer]

@shreyas02 I already did this for non-transient files in PR #158 . I think the result is quite more elegant. I would say it's just a matter to applying the same treatment to the necessary functions.

Note this checks that the triangulations are there but empty. If you are working within a sub-communicator, you should just not be trying to print from ranks that are empty. You can easily create a new subcomm using the generate_subparts functionality. You can have a look at how I deal with this in GridapSolvers.

[shreyas02]

Thanks @JordiManyer! However, there is MPI communication across MPI processes taking place multiple times in your implementation. Is this okay?

https://github.com/gridap/GridapDistributed.jl/blob/3e4fafec239ba98000839f6affb8f0dcb69217c3/src/Visualization.jl#L184C1-L192C4

Here is the implementation using this function -

struct DistributedPvd{T<:AbstractArray}
  pvds::T
  parts::AbstractArray
end

function Visualization.createpvd(trian::DistributedTriangulation,parts::AbstractArray,args...;kwargs...)
  nparts, new_parts = filter_empty_parts(parts,local_views(trian))
  pvds = map(new_parts) do part
    if part == 1
      paraview_collection(args...;kwargs...)
    end
  end
  DistributedPvd(pvds,new_parts)
end

function Visualization.createpvd(f,trian::DistributedTriangulation,parts::AbstractArray,args...;kwargs...)
  pvd = createpvd(trian,parts,args...;kwargs...)
  try
    f(pvd)
  finally
    savepvd(pvd)
  end
end

function Visualization.createpvd(parts::AbstractArray,args...;kwargs...)
  pvds = map(parts) do part
    if part == 1
      paraview_collection(args...;kwargs...)
    end
  end
  DistributedPvd(pvds,parts)
end

function Visualization.createpvd(f,parts::AbstractArray,args...;kwargs...)
  pvd = createpvd(parts,args...;kwargs...)
  try
    f(pvd)
  finally
    savepvd(pvd)
  end
end

function Visualization.savepvd(pvd::DistributedPvd)
  map(pvd.pvds, pvd.parts) do pvd, part
    if part == 1 
      vtk_save(pvd)
    end
  end
end

function Base.setindex!(pvd::DistributedPvd,pvtk::AbstractArray,time::Real)
  map(pvtk) do pvtk
    if !isnothing(pvtk)
      vtk_save(pvtk)
    end  
  end
  map(pvtk,pvd.pvds,pvd.parts) do pvtk,pvd,part
    if part == 1
      pvd[time] = pvtk
    end
  end
end

[shreyas02]

Also, @JordiManyer, PR #169 breaks my code.

[JordiManyer]

However, there is MPI communication across MPI processes taking place multiple times in your implementation. Is this okay?

The advantage is that you get contiguous numbering of the files (which I think is necessary? I dont completely remember....).
This is only activated when saving stuff, I think we do not care about performance that much.

[JordiManyer]

Also, @JordiManyer, PR #169 breaks my code.

This is so unspecific. Unless you give me some context, I cant really help you....

Considering most of what I touched is completely new, I assume the issue is in the get_cell_gids changes?

[shreyas02]

@JordiManyer sounds good. Do I create a PR with the proposed changes for handling the pvd files.

I’ll probably start a new issue regarding the other issue. (#172)

[shreyas02]

Hi @JordiManyer, shall I add the draft that I sent above handling pvd files as a pr?

[JordiManyer]

Sure. Is the struct new? If so, is it necessary?

[shreyas02]

Yes, it is new. I need the parts abstract array later. For example, this is how I use it -

map(pvd.pvds, pvd.parts) do pvd, part
    if part == 1 
      vtk_save(pvd)
    end
  end

Also, there are two methods for creating a pvd now -

1. createpvd(parts::AbstractArray,args...;kwargs...) (old )
2. createpvd(trian::DistributedTriangulation,parts::AbstractArray,args...;kwargs...) (newly added)