How to read all the dof elements in a cell in Petsc MPI problem

[Rahul Gopalan Ramachandran]

 Hello everyone,

I am having problem with applying static condensation as explained in step-44 in parallel. The objective is to employ a parallel solver. Any guidance on the issue would be much helpful.

As of now, the code works when run on single MPI process, but fails when np > 1. The error originates from the following part of the code where it tries to read the elements of the tangent_matrix, which is a PETScWrapper::MPI::BlockSparseMatrix. The issue arises trying to access the dofs not owned by the MPI process at the boundary partition. 

//-----------------------------------------------------------------------------

for (unsigned int i=0; i<dofs_per_cell; ++i)
      for (unsigned int j=0; j<dofs_per_cell; ++j)      
        data.k_orig(i,j) =  tangent_matrix.el(data.local_dof_indices[i],
                            data.local_dof_indices[j]); 

//-----------------------------------------------------------------------------

In the case of  PETScWrapper::MPI::BlockVector I can circumvent such a problem defining a locally_relevant BlockVector as below. I am able to print out all the elements in this way (I am not sure if this is the most efficient way!).

//-----------------------------------------------------------------------------

// Access ghost elements in system_rhs
    LA::MPI::BlockVector system_rhs_local (locally_owned_partitioning,
                                         locally_relevant_partitioning,
                                         mpi_communicator);
    system_rhs_local = system_rhs;
    // printing out all the system_rhs values for all dofs in the cell
    for (unsigned int i=0; i<dofs_per_cell; ++i)
        if (this_mpi_process == 1)
           std::cout << " \n Cell:"
              << cell
              << ";system_rhs("<< i <<"):"
              <<system_rhs_local(data.local_dof_indices[i])
              << " ; this_mpi_process: "
              << this_mpi_process << std::flush;

//-----------------------------------------------------------------------------

Below is the initialization of the tangent_matrix BlockSparseMatrix for reference.

//-----------------------------------------------------------------------------

BlockDynamicSparsityPattern dsp(locally_relevant_partitioning);
DoFTools::make_sparsity_pattern(dof_handler, coupling, dsp, constraints, false);
  

SparsityTools::distribute_sparsity_pattern(
    dsp,
    dof_handler.locally_owned_dofs(),
    mpi_communicator,
    locally_relevant_dofs);

    tangent_matrix.reinit(locally_owned_partitioning,dsp, mpi_communicator);

//-----------------------------------------------------------------------------

Please let me know if any more details are needed.

Looking forward to your input.

Thank you and Regards,

Rahul G R

Post Doc - MPI PKS, Dresden

[Wolfgang Bangerth]

On 9/27/22 11:54, Rahul Gopalan Ramachandran wrote:
>
> I am having problem with applying static condensation as explained in step-44
> in parallel. The objective is to employ a parallel solver. Any guidance on the
> issue would be much helpful.
>
> As of now, the code works when run on single MPI process, but fails when np >
> 1. The error originates from the following part of the code where it tries to
> read the elements of the tangent_matrix, which is a
> PETScWrapper::MPI::BlockSparseMatrix. The issue arises trying to access the
> dofs not owned by the MPI process at the boundary partition.
>
> //-----------------------------------------------------------------------------
> for (unsigned int i=0; i<dofs_per_cell; ++i)
>       for (unsigned int j=0; j<dofs_per_cell; ++j)
>         data.k_orig(i,j) =  tangent_matrix.el(data.local_dof_indices[i],
>                             data.local_dof_indices[j]);
> //-----------------------------------------------------------------------------
>

Rahul -- the problem is that, unlike vectors, matrices are generally not
written to use data structures that provide you with the equivalent of "ghost
elements". That is because, in most applications, matrices are "write only":
you build the matrix element-by-element, but you never read from it
element-by-element. As a consequence, people don't provide the facilities to
let you read elements stored on another process.

If you need this kind of functionality, you probably have to build it
yourself. The way I would approach this is by looking at your algorithm,
figuring out which matrix elements you need to be able to read for rows not
locally owned (most likely the rows that correspond to locally-relevant but
not locally-owned; it may also be locally-active but not locally-owned). You
will then have to read these elements on the process that owns them and send
them to all processes that need them but don't own them. In the current
development version, you can use functions such as Utilities::MPI::Isend to do
that with whatever data structure you find convenient, but it can also be done
with standard MPI calls.

This isn't particularly convenient, but it is the best anyone can offer short
of writing the necessary functionality in PETSc itself.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Rahul Gopalan Ramachandran]

Hello Dr. Bangerth,

Thank you for the clarification. Isend should be the way to patch this. Alternatively, what is your opinion on copying the matrix to a single processor (if that is possible) using MPI_Gather? Would the cell iterators then no longer work? Ignoring scalability will this also be an option?

Regards,
Rahul

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> --- You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/3PXbitrTj_k/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/6c1a9f21-2918-37ac-87e8-fd0fe8614bdb%40colostate.edu.

[Wolfgang Bangerth]

On 9/27/22 23:50, Rahul Gopalan Ramachandran wrote:
> Thank you for the clarification. Isend should be the way to patch this. Alternatively, what is your opinion on copying the matrix to a single processor (if that is possible) using MPI_Gather? Would the cell iterators then no longer work? Ignoring scalability will this also be an option?

At least in principle, of course, we'd like to avoid writing programs
that we know can't scale because each process stores data replicated
everywhere -- like the entire matrix. In practice, if your goal is to
run on 10 or 20 processes, this may still work, though you should
recognize that the system matrix is the largest object you probably have
in your program (even if you fully distribute it).

If your goal is to replicate the matrix on every process, it might be
easiest to create a data structure that collects all of the local
contributions (local matrix, plus the dof indices array) and send that
around between processes, which then all build their own matrix. That
may be simpler than sending around the matrix because the latter is a
PETSc object to which you don't easily have access to its internal
representation.

[Rahul Gopalan Ramachandran]

> "If your goal is to replicate the matrix on every process, it might be easiest to create a data structure that collects all of the local contributions (local matrix, plus the dof indices array) and send that around between processes, which then all build their own matrix. That may be simpler than sending around the matrix because the latter is a PETSc object to which you don't easily have access to its internal representation."

Thank you Dr. Bangerth for this suggestion. As I am fairly new to object oriented programing, I did not think in this light. I have a related question.May I ask if there is a function in dealii to identify cells on the boundaries of subdomains?

Regards,
Rahul

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> --- You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/3PXbitrTj_k/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/eb258ba9-a221-5229-17e5-61453b639cf6%40colostate.edu.

[Wolfgang Bangerth]

On 9/29/22 01:21, Rahul Gopalan Ramachandran wrote:
> Thank you Dr. Bangerth for this suggestion. As I am fairly new to object oriented programing, I did not think in this light. I have a related question.May I ask if there is a function in dealii to identify cells on the boundaries of subdomains?

There are a number of functions in namespace GridTools that can help
with this. Specifically, the following come to mind:
GridTools::compute_ghost_cell_halo_layer
GridTools::exchange_cell_data_to_ghosts
GridTools::compute_vertices_with_ghost_neighbors