Dealii: Refactor interfaces to parallel vectors in MatrixFree/FEEvaluation

[davydden]

1. Currently matrix_free.h contains some internal functions update_ghost_values_start(), reset_ghost_values(), update_ghost_values_finish(), compress_start(), compress_finish() which need to be re-implemented for each new vector type to be used with matrix-free. Should not we move those to VectorSpaceVector or alike? It would probably also make sense to move MPI partitioner to VectorSpaceVector. Then a lot of things will be easier.

2. Similarly, i would deprecate MatrixFree<>::initialize_dof_vector() as it also accumulates custom definitions inside matrix_free.h for each new vector class. I think we just need to ask to use MPI partitioner from MatrixFree data object to initialize anything derived from VectorSpaceVector, i.e. LinearAlgebra::distributed::Vector can be initialized using MPI paritioner. That would also make it more clear what interface vectors need to provide/implement in order to be used with MatrixFree.

3. Similar things happen in fe_evaluation.h with internal vector_access() functions which give read/write access to an element of a vector and check_vector_compatibility(). Currently it also necessitates adding extra code for each new vector to be used with FEEValuation. Those probably should be related to ReadWriteVector.

I know that @kronbichler has some change to Matrix-Free in the pipeline and probably most of those issues are addressed, but I nonetheless wanted to document somewhere what I perceive as limitations, or rather inconveniences.

[kronbichler]

Let me add a few remarks here:

Should not we move reset_ghost_values(), update_ghost_values_finish(), compress_start(), compress_finish() from matrix_free.h to VectorSpaceVector or alike?

While I agree that we need to find a way to express the update in the right way, I do not think that we should move the update ghost values to the vector space vector. The idea is to created ghosted vectors as ReadWriteVector classes. As discussed in #2535, I would like to be able to create lightweight objects that do not copy any data but simply attach to the data array in LinearAlgebra::distributed::Vector.

[kronbichler]

Sorry, I pressed the wrong button... Here my post continues:
Given the different states, a ReadWriteVector class would be the natural use inside MatrixFree::loop contexts. Alternatively, we might completely hide the vectors once they are fed into the loop. How could we do that? I'm thinking of attaching one (or more) FEEvaluation-type object to each vector. This would avoid having to pass in block vectors or an std::vector of vectors as well. But this means that we need to re-work the way the evaluators are set up right now, which is necessary to some extent with the GPU extension developed by @kalj and @Rombur anyway, so we should try to get it right once.

On the other hand, I do see the problem that we cannot completely hand over the update task to the vectors, at least not exclusively. MPI communication is not for free and one would like to keep it down to some basic amount. In DG applications, you might apply mass-matrix type operations (no communication necessary) on the same vector as when you compute face integrals that cross processor boundaries. A have ideas of how to express that, but in order to organize the data exchange I need to get that from the vector to MatrixFree. We could possibly do that just for one vector type LinearAlgebra::distributed::Vector because it needs some particular steps. I'm about writing these things down and will let you have a read of the why and how.

Similarly, i would deprecate MatrixFree<>::initialize_dof_vector() as it also accumulates custom definitions inside matrix_free.h for each new vector class. I think we just need to ask to use MPI partitioner from MatrixFree data object to initialize anything derived from VectorSpaceVector, i.e. LinearAlgebra::distributed::Vector can be initialized using MPI partitioner.

I definitely agree on that. We initially didn't have the vector partitioner organization completely clear but I agree that we should do it that way. It simplifies things.

Similar things happen in fe_evaluation.h with internal vector_access() functions

As you say, this is linked with the way the MPI update works. Right now, it goes through some specializations and the bad news is that DG-type codes and some optimizations make things probably worse because you want to do vectorized read operations where you basically only want a pointer to the vector's data array vector.begin(). I know that your work necessitates dedicated local_element(i) functions, so we would need two interfaces at least, but probably not more than that. When we re-work these things, we should also make sure to also fail on vectors that do not support this model - I think we currently fail somewhere deep inside the access functions when we have distributed Trilinos/PETSc vectors.

[kronbichler]

Another thing we need to improve is the handling of block vectors and components of systems (this is related to #2000): Right now we only support FESystem consisting of a single element type, but we really want to come to the point where we also support multiple elements, say Taylor-Hood, within a single DoFHandler. As far as the evaluators are concerned, this is simple because one simply needs to add several internal::MatrixFreeFunction::DoFInfo objects, one for each base elements. One thing that is more complicated is that the definition of a BlockVector in the MPI case does not currently go together with LinearAlgebra::distributed::Vector and the related Utilities::MPI::Partitioner classes that only support a single local range. I think this should be done by appropriate storage in DoFHandler.

[kronbichler]

Even though it would have been nice to make a step forward on this topic, the time constraints are such that we will only be able to fix #5667 for the 9.0 release. Adjusting the target milestone.

[kronbichler]

Partly addressed by #7782; moving target milestone regarding a more general solution.