Extracting a sub-block matrix (vector-valued problem)

[Andrey Volkov]

Hello everyone,

Let us consider a vector-valued problem, which weak form of the discrete equations results in the following linear system:

While this system could be solved as a whole, let us consider the case when we partition the global block matrix into "sub-block matrices" as follows:

Now, we can first solve the first sub-block matrix (one can use the Schur complement): 

and then solve the second one:

So, how to extract a submatrix from the global block matrix? In my view, there are two approaches:

A) Let us assume we have already constructed the global matrix, and its blocks have been established. The global matrix has a sparsity pattern represented by BlockDSP of the type BlockDynamicSparsityPattern. Now we want to initialize the sub-block matrices, but to do that, we need to have the sparsity patterns for each of the "sub-problems". Let us initialize the "sub-sparsity" pattern for the first sub-block and add non-zeros entries to it based on the global sparsity pattern as follows:

//--- creating the BDSP for the fluid flow problem

BlockDynamicSparsityPattern BlockDSP_fluid_flow(dofs_per_block_fluid_flow, dofs_per_block_fluid_flow);

//--- copying the block of the sparsity pattern that rperesents only fluid flow problem

for(unsigned int i = 0; i < total_dofs_fluid_flow; i++)

    for(unsigned int j = 0; j < total_dofs_fluid_flow; j++)

        if (this->BlockDSP.exists(i,j))

            BlockDSP_fluid_flow.add(i,j);

//--- converting BDSP to BSP

BlockSparsityPattern sparsities_fluid_flow;

sparsities_fluid_flow.copy_from(BlockDSP_fluid_flow);

where the  dofs_per_block_fluid_flow is a vector of DoFs per each block of the first sub-problem (sizes of the blocks M and B),  total_dofs_fluid_flow is the total number of DoFs of the first sub-problem.

Thus, we have extracted the sparsity pattern for the first sub-problem, and we can now create the corresponding sub-block matrix:

//--- initializing the matrix block for fluid flow problem only

BlockSparseMatrix<double> matrix_fluid_flow(sparsities_fluid_flow);

Finally, we fill  matrix_fluid_flow with the corresponding values from the global block matrix as follows:

//--- copying content of the total matirx to the matrix that is specific to the fluid flow problem

for(unsigned int i = 0; i < total_dofs_fluid_flow; i++)

    for(unsigned int j = 0; j < total_dofs_fluid_flow; j++)

        if ( sparsities_fluid_flow.exists(i,j) )

            matrix_fluid_flow.set(i, j, this->matrix.el(i,j) );

where matrix is a global block matrix.

Then we use linear solvers to obtain a solution for the first sub-problem, and after that, we do the same extraction/solution procedure for the second sub-problem. Thus, we have solved the two sub-systems separately.

I am sure at this point you notice that iterating through all the entries of the global matrix and its sparsity pattern is quite an expensive operation. It seems like solving sub-systems separately should, in fact, improve the stability and speed, but going through all these for-loops increased the simulation time.

B) The second approach is to never construct the global block matrix and instead create two block matrices based on their sparsity patterns and fill it with the values using separate routines. Unfortunately, this would involve rewriting a lot of the code that I did not write (I am working on the OpenFCST).

So, before doing this, I would like to inquire the following questions:

• Is there a functional in the deal.ii that allows extracting mentioned above "sub-block matrix" from the global block matrix? Could recommend something?
• Or perhaps I should go on and try to implement something similar to the .block(i,j) method of the BlockMatrix class, but that allows to access the multiple blocks (somehow based on the sub_object of the BlockMatrix)?

Hopefully, the above explains my problem well, but please let me know if you need further clarification. I will appreciate any help, and thank you for creating a great FEA library and constantly improving it!

If you are curious, I am trying to solve the linearized formulation (i.e. Newton-Raphson method) of the incompressible Navier-Stokes problem coupled with multiple reaction-diffusion equations.

---

Respectfully,

Andrey Volkov

[Wolfgang Bangerth]

On 3/14/21 2:33 PM, Andrey Volkov wrote:
>
> So, before doing this, I would like to inquire the following questions:
>

> * Is there a functional in the deal.ii that allows extracting mentioned

> above "sub-block matrix" from the global block matrix? Could recommend
> something?

Yes:
system_matrix.block(0,0)
gives you a reference to the top left block :-)

You probably want to take a look at how step-20 or step-22 solve their (block)
linear systems. They use exactly the kind of approach you are looking for.

Best
W.


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

[Andrey Volkov]

Dear Dr. Bangerth,

The system_matrix.block(0,0) gives access to a single block, namely M (correct?). However, is there a way to access/extract multiple blocks as a whole block matrix, let us say { [M B], [B 0] } or { [D E], [F G] }?

Sure, I looked at step-20/22, as all I have done so far is based on these tutorials.

Thanks for a quick reply.

[Wolfgang Bangerth]

On 3/14/21 3:23 PM, Andrey Volkov wrote:
>
> The system_matrix.block(0,0) gives access to a single block, namely *M
> *(correct?). However, is there a way to access/extract _multiple blocks as a
> whole block matrix_, let us say { [M B], [B 0] } or { [D E], [F G] }?

>
> Sure, I looked at step-20/22, as all I have done so far is based on these
> tutorials.

Not quite so easily, but ASPECT does this: It simply defines its own 2x2 block
matrix class that points to sub-blocks of a larger matrix. See here:

https://github.com/geodynamics/aspect/blob/master/source/simulator/solver.cc#L47-L105

Of course, if you don't actually need access to the M and B blocks
individually, you don't have to structure your matrix that way. The
structuring of the matrix based on solution variables is your choice, and you
can have a 2x2 block pattern for a matrix even if you have 3 solution
variables, and then you can put [[M, B],[B, 0]] all into one block.

[Andrey Volkov]

Well, yes, I'd like to have access to individual blocks as well (i.e., M and B), but anyway, I know what to do now.

Thanks for the example and helpful discussion.