Issues found/lessons learned with MPI-parallel Stokes solver

[amartinhuertas]

In this issue, I summarize the issues/lessons that I found/learned while developing MPIPETScDistributedStokesTests.jl, a parallel MPI program for the solution of the Stokes PDE. In order to avoid building a pressure FE space composed of zero mean functions (for reasons made clear below), the script, as of today (see link above), is solving a problem with manufactured solution on a square, with Dirichlet boundary conditions everywhere, except for the lower edge of the square, where we have Neumann boundary conditions. The issues are (prioritized by importance):

• [UPDATE: Solved in 2a71fe7 1663849 c754797 0b95e0f 94df2cd e92b9c3] We are currently using a software design pattern in which we have a DistributedFESpace composed of local FE spaces. The polymorphism is in the local FE spaces (e.g., SingleFieldFESpace versus MultiFieldFESpace). DistributedFESpace is unique. DistributedFESpace glues together the parts by building a global DoF numbering across processors. While this pattern has quite a lot of expressivity, I do not see how it can be used to cleanly build a global zero mean FE space. We need some sort of global layer of code in which the processors agree in the single DOF to be fixed, and the zero mean constraint imposition has to be computed globally by computing the pressure mean/domain volume by assembling (reduce-sum) local contributions from each of the subdomains. Is it possible to acommodate in the interplay of DistributedFESpace+Local FE spaces these operations? Do you foresee any other similar scenario? (e.g., FESpaces with multipoint linear constraints?)

• [UPDATE: Solved in 2a71fe7 1663849 c754797 0b95e0f 94df2cd e92b9c3] Just guessing, the problem in the point above could be solved if we had at least two types of SingleField DistributedFESpaces (the standard one, and the one that imposes the constraint), and a MultiFieldFEspace defined as the composition of global FE spaces, but this clearly breaks our original intended design pattern of having a single DistributedFEspace and besides, it needs extra communication (see issue Misc tasks (in the lack of a better name) #3).

• I had to define a new num_dirichlet_dofs() method for MultiFieldFESpace, even this FESpace is not a SingleFESpace. This is is required in this line.

• Also worths mentioning a very annoying BUG that drove me crazy during many hours, and that I consider convenient to document here in case it is useful for others. In particular, when restricting the multifield trial basis into its velocity and pressure components within the function that defines the bilinear form (see this line, if I use u as the variable name for velocity (note that the code is currently using y), then this variable assignment has a side-effect in the local scope of the run function, in particular, it modifies the binding of the local variable u of run previously set up in this line. This causes this line to crash later on as u is no longer an object of type Function.

• [UPDATE: This BUG has been already solved by @fverdugo in Gridap.jl] We need a bug-fix for the following issue in Gridap.jl: AssertionError triggered in the call path rooted at AppliedArrays: BUG? Gridap.jl#323

• [UPDATE: This issue has been solved after the excellent work by @fverdugo in PR #376.] In order to let the method that generates the global DoF numbering (see here) work with local MultiFieldFESpaces, I had to modify the current implementation of MultiFieldArray{T,N} in Gridap.jl such that it now extends AbstractArray{T,N}. The particular lines in which this is required are this one and this one. This has been implemented in a branch of Gridap.jl. See here for the current implementation of MultiFieldArray{T,N}. However, following concerns/issues are around:

1. For obvious reasons, getindex() semantics were changed from block-based semantics to scalar-based semantics, i.e., A[1,1] does not refer to the upper left block in the block structure of A, but the entry in the upper left corner of A. This already causes some Gridap.jl tests to fail. See here for details.

2. I added two member variables in order to accelerate indexing operations of MultiFieldArray{T,N}. However, there are case uses of this data type in which the blocks member variable is set up with undef (see here), and thus these member variables cannot be set up during creation, and cannot be created later on as the struct is immutable. In conclusion, I am not quite happy with the current solution.

3. Again guessing, I have the impresison that MultiFieldArray{T,N} is a multi-purpose data structure, that is re-used in many scenarios. A blocked cell dofs array is just one of these scenarios. Perhaps we could create a separate ad-hoc data type for such purpose. This would solve 1. but 2. would still be there. Any ideas?