Initialization of PETSc::MPI::BlockSparseMatrix

[Gabriel Stankiewicz]

Dear deal.II  Users,

I am trying to setup an eigenvalue problem of a piezoelectric bimorph cantilever, in which I need to use a distributed block sparse matrix. My problem contains displacement and electric potential degrees of freedom. Because of that, I get 0 entries in the mass matrix that correspond to the potential DoFs. My approach is to condense the stiffness and mass matrices so that I can get rid of the potential DoFs in the system. For that, I decided to use PETSc::MPI::BlockSparseMatrix.

When I try to initialize the PETSc::MPI::BlockSparseMatrix, I pass an instance of std::vector<IndexSet> that contains the locally owned DoFs for the displacement block and the electric potential block. However, the reinit() function:

void BlockSparseMatrix< number >::reinit
(
const std::vector< IndexSet > & 
rows,
const std::vector< IndexSet > & 
cols,
const BlockDynamicSparsityPattern & 
bdsp,
const MPI_Comm & 
com  )

fails to pass the assertion that compares the individual block sizes and the sizes of the IndexSets. Since the size of an IndexSet needs to correspond to the total number of DoFs in the system (otherwise I will not be able to store some of the DoF indices), but the block sizes are naturally smaller than the total number of DoFs, this assertion cannot be passed.

The violated condition was:
    rows[r].size() == bdsp.block(r, c).n_rows()
Additional information:
    invalid size

Could anyone tell me if my understanding is correct or if there is a solution to my problem?

Best regards,

Gabriel Stankiewicz

[Timo Heister]

Hi,

I don't quite understand your reasoning for a block with 0 DoFs. Try removing any rows/columns from your block sparsity pattern and from your IndexSet's that have size 0. Notice that the number of blocks is independent of the number of components in your FEM.

Can you post the values for rows[r].size() and bdsp.block(r,c).n_rows() for all rows r? 

Best,

Timo

[Gabriel Stankiewicz]

Hi Timo,

thanks for replying. I think I didn't describe my problem correctly. The 4 blocks have sizes: 1815 x 1815, 1815 x 330, 330 x 1815, 330 x 330. Such a division results from having 1815 displacement DoFs and 330 electric potential DoFs in the eigenvalue problem. Only the block 1815 x 1815 of the mass matrix has non-zero entries and there are 330 zeros in the diagonal of the mass matrix. I would like to condense the system to get rid of the zero diagonal entries in the mass matrix using the operations on the blocks. 

Since the total number of the DoFs in the system is 2145, all the IndexSets, that are passed to the reinit() function of the PETSc::MPI::BlockSparseMatrix are of size 2145. 

This results in failed assertion, since: rows[0].size() = 2145, rows[1].size() = 2145, but bdsp.block(0,0).n_rows() = 1815, bdsp.block(1,1).n_rows() = 330 etc...

Sorry, I cannot post the screenshots of these values since I am on holidays starting from today, I remember them from previous days as I checked them already at that time.

I tried to explicitly set the size of the IndexSets, i.e. for the IndexSet corresponding to electric potential DoFs that would be 330. Then, I cannot add indices higher than 330 to this IndexSet due to its internal structure.

Is there a way to overcome that or is the assertion itself wrong?

Best regards,

Gabriel

[Timo Heister]

> This results in failed assertion, since: rows[0].size() = 2145, rows[1].size() = 2145, but bdsp.block(0,0).n_rows() = 1815, bdsp.block(1,1).n_rows() = 330 etc...

The size() of each IndexSet needs to correspond to the size of each
block of the sparsity pattern and this is what the assertion checks. I
think that this check is correct. If you want your first block to be
1815 DoFs large, then rows[0].size() needs to be 1815.

How are you creating these IndexSets? I find the function
IndexSet::split_by_block() convenient to use.


--
Timo Heister
http://www.math.clemson.edu/~heister/

[Gabriel Stankiewicz]

I created IndexSets by using the function DoFTools::locally_owned_dofs_per_component() and then gathering all indices corresponding to displacement DoFs (three instances of IndexSet) into one IndexSet using IndexSet::add_indices() and the fourth instance correponded to electric potential DoFs. However, the function DoFTools::locally_owned_dofs_per_component() wasn't returning only locally owned DoFs but also the locally not owned ones. In the end I used the check locally_owned_dofs.is_element() to make sure I gather only locally owned DoFs for displacement and electric potential into my final IndexSets. The resulting IndexSets were always of size 2145. 

I also tried first explicitly setting the sizes of the final IndexSets using IndexSet::set_size() function, but as said, I could not add indices larger than the size to these index sets (eg. id = 440 could not be added to the IndexSet of size 330).

I didn't try the function IndexSet::split_by_block(), I didn't notice it in the dealii manual. I will try this one when I come back to the office, thank you for the suggestion!

Best,

Gabriel

[Wolfgang Bangerth]

On 9/17/20 9:07 AM, Gabriel Stankiewicz wrote:
> I created IndexSets by using the function
> DoFTools::locally_owned_dofs_per_component() and then gathering all indices
> corresponding to displacement DoFs (three instances of IndexSet) into one
> IndexSet using IndexSet::add_indices() and the fourth instance correponded to
> electric potential DoFs. However, the function
> DoFTools::locally_owned_dofs_per_component() wasn't returning only locally
> owned DoFs but also the locally not owned ones.

No, that isn't right. This is a function that is so widely used that it is
hard to believe that it would be wrong. If you think that it returns something
wrong, I believe that your *understanding* of what this function returns is
wrong, and that might lead to further confusion downstream when you are
thinking how to initialize matrices.

There are quite a number of tutorial programs that build parallel block
matrices. I know that step-32 does this for Trilinos, but I think that the
parallel Navier-Stokes solver does it for PETSc. It would be useful to just
compare your code with how these programs do it.

Best
W.

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

[Gabriel Stankiewicz]

Thank you for your input!

1. Regarding the initialization of PETSc::MPI::BlockSparseMatrix:

I have used the the IndexSet::split_by_block() function and this indeed works good. Thanks for the suggestion!

Unfortunately, I have encountered another issue. The PETSc::MPI::BlockSparseMatrix must be partitioned into contiguous chunks, for that reason renumbering the DoFs by component in order to define blocks for displacement DoFs and electric potentials DoFs fails, giving an error: "PETSc only supports contiguous row/column ranges". I know that Trilinos allows for non-contiguous partitioning (according to what is written in step-32), but I need PETSc to run the eigenvalue problem using SLEPc. Do you have any ideas how this issue could handled?

2. Clarification to locally_owned_dofs_per_component(), which I mentioned in my previous post:

I created a small working example to demontrate what I described in the previous post (see attached). I run the file using "mpirun -np 2 test_dofs" in debug mode. Please check the attached image for the output.

This shows that the function locally_owned_dofs_per_component() divides the DoFs correctly per component but not per processor. According to what is written in the documentation, the union of locally owned DoFs for all components should correspond to locally owned DoFs, which is not the case in this example.

Best regards,

Gabriel

[Wolfgang Bangerth]

Gabriel,

> 1. Regarding the initialization of PETSc::MPI::BlockSparseMatrix:
>
> I have used the the IndexSet::split_by_block() function and this indeed works
> good. Thanks for the suggestion!
>
> Unfortunately, I have encountered another issue. The
> PETSc::MPI::BlockSparseMatrix must be partitioned into contiguous chunks, for
> that reason renumbering the DoFs by component in order to define blocks for
> displacement DoFs and electric potentials DoFs fails, giving an error: "PETSc
> only supports contiguous row/column ranges". I know that Trilinos allows for
> non-contiguous partitioning (according to what is written in step-32), but I
> need PETSc to run the eigenvalue problem using SLEPc. Do you have any ideas
> how this issue could handled?

We are all using PETSc all the time, so that can't be true.

In the end, you need a block matrix each block of which is partitioned between
processors. There are numerous tutorial programs that do exactly this. Why
don't you copy the way step-32, for example, enumerates and partitions degrees
of freedom? There, all that is done is to call

stokes_dof_handler.distribute_dofs(stokes_fe);

std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
stokes_sub_blocks[dim] = 1;
DoFRenumbering::component_wise(stokes_dof_handler, stokes_sub_blocks);

and that is good enough to have things partitioned into blocks, and have
everything within each block be contiguous.

> 2. Clarification to locally_owned_dofs_per_component(), which I mentioned in
> my previous post:
>
> I created a small working example to demontrate what I described in the
> previous post (see attached). I run the file using "mpirun -np 2 test_dofs" in
> debug mode. Please check the attached image for the output.
>
> This shows that the function locally_owned_dofs_per_component() divides the
> DoFs correctly per component but not per processor. According to what is
> written in the documentation, the union of locally owned DoFs for all
> components should correspond to locally owned DoFs, which is not the case in
> this example.

The problem here is that you are using a sequential triangulation. On every
processor, the DoFHandler owns all DoFs and that's what you see when you
output the DoFs per component: every dof is owned by every process, because
they simply don't know about each other. (The fact that you partition, i.e.
set subdomain_ids, on the triangulation doesn't matter: You are using the
sequential triangulation class, it knows nothing of other processors.)

On the other hand, you call DoFTools::locally_owned_dofs_per_subdomain which
simply counts DoFs per subdomain -- but for a sequential triangulation, the
subdomain_id of every cell is just a number without any special meaning.

[Gabriel Stankiewicz]

Dear Prof. Bangerth,

Regarding the locally_owned_dofs_per_component() function:

Thank you for the detailed explanation, this clarifies the problem. Now I have a better understanding of how a sequential triangulation is handled when partitioned and used for problems ran on multiple processors. The confusion came from the fact that printing the elements of locally_owned_dofs actually did output only the "locally owned" ones, but the locally_owned_dofs_per_component output "all" DoFs.

Regarding the PETScWrappers::MPI::BlockSparseMatrix:

I actually did setup the matrix following the step-32, but found that the error that the partitions are not contiguous is also related to the sequential triangulation. As you described, the DoFHandler considers every cell as locally owned, as opposed to parallel::shared::Triangulation to which I switched now and do not get the error message. See the attached example, which works for me now.

Thank you again both Timo Heister and Prof. Bangerth for your time and help.

Best regards,

Gabriel