New block array implementation

[santiagobadia]

After talking to @amartinhuertas, @fverdugo I think that we probably need to define a new implementation of BlockArray.

The current BlockArray is an array of arrays with some metadata (touch). I think it is ok for same dimension multiphysics (if you do not want to statically condense interior DOFs). However, when moving to mixed dimensions, I think this is not the best way to go. In hybrid formulations, it is mandatory to locally condense interior DOFs. Interior DOFs can involve multiphysics. As a result, it would be great if the interior DOFs block would be a dense matrix. We could use level-3 BLAS to optimise this (computationally very intensive) step. E.g., if 1 denotes interior and 2 denotes interface, we need to compute A22 - A21*inv(A11)*A12. With the current approach, this is not possible.

One sub-optimal solution is to copy the BlockArray A11 into a new dense matrix. But I think this is inefficient. We would also need to copy A12 and A21 to dense matrices to use level-3 BLAS. We could pursue this line as an initial solution and try to go further. Once we will have hybrid methods working, we know what we can improve.

A final solution (IMO) is to define an AbstractBlockArray class, with the current struct as a sub-type, e.g., BlockArraySparse, and another struct BlockArrayDense (working with a dense matrix). We just need to implement the current queries working with offset. It does not seem :rocket: science. The idea is that we can also combine them recursively. we should also need to know how to efficiently implement a broadcast between these new blocked matrices.

What do you think?

[amartinhuertas]

For completeness, let me attach some hand-written (quick & dirty) notes that I wrote during the meeting to illustrate the situation. I am going to think a little bit more about the proposal in the previous post to see whether I can identify additional pros & cons (if any). BTW, the data type at hand is called ArrayBlock not actually BlockArray (I know that it is not a big deal, but I guess that BlockArray is already a data type offered by an existing package https://github.com/JuliaArrays/BlockArrays.jl).

2021-05-21-Note-11-38.pdf .

[amartinhuertas]

Let me add the following remark to @santiagobadia's post (while I am still thinking). I think that the current implementation of ArrayBlock is well-suited (best-suited?) for blocked arrays of Fields. It is in the evaluation of the current ArrayBlock on a set of points where the limitation is exposed (as we are using the same BlockMap/ArrayBlock pair for such purpose). Thus, the solution could come using (say) ArrayBlockDense and MapBlockDense only when evaluating ArrayBlocks of Fields on a set of points and performing a broadcasted operation among them. Still I have to confirm/understand that using ArrayBlockDense we would not end up with temporary block matrices plenty of zeros. I am positive that this would not happen, due to the way LazyArrays are combined, etc.), but I am on the way to confirm this by digging deeper into the code.

[amartinhuertas]

To complement my previous post with examples:

• The following line creates a blocked array of Fields out an existing array of Fields. https://github.com/gridap/Gridap.jl/blob/9b3c7fbe63730c09cfed052eebc62afb9bf5de2c/src/MultiField/MultiFieldFESpaces.jl#L124
• The following line "evaluates" a blocked array of Fieldss on a set of points: https://github.com/gridap/Gridap.jl/blob/9b3c7fbe63730c09cfed052eebc62afb9bf5de2c/src/Fields/ArrayBlocks.jl#L290

[fverdugo]

E.g., if 1 denotes interior and 2 denotes interface, we need to compute A22 - A21inv(A11)A12.

@santiagobadia can you write the static condensation expression also involving field ids? It is not clear for me if say A22 contains DOFs of all fields or only one.

[fverdugo]

@amartinhuertas @santiagobadia when you refer to "dense" matrices, do you mean to add rows and cols of the different fields in the same dense matrix? If so, it is likely that this would re-introduce the problems I fixed in the last refactoring where I implemented ArrayBlock. Apart that you will end up with matrices full of zeros.

So I would like to understand if it really needed to use dense matrices.

[santiagobadia]

All fields. We need to compute the inverse of the multiphysics block, otherwise it would be wrong.