MGLevelGlobalTransfer/MGTransferMatrixFree: allow user to init vectors

[peterrum]

Similar to:

dealii/include/deal.II/multigrid/mg_transfer_global_coarsening.h

Lines 353 to 356 in 15440af

	MGTransferGlobalCoarsening(
	const MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> &transfer,
	const std::function<void(const unsigned int, VectorType &)>
	&initialize_dof_vector = {});

Not completely finished. I'll continue tomorrow. Hopefully, this is the start of the end of adjust_ghost_range_if_necessary() and the derived MGTransferMatrixFree operators that only extend copy_to_mg().

[kronbichler]

I have not looked into the code yet, but this is a good help. We should keep an eye on whether or not this can be combined with the code currently used for a similar purpose

dealii/include/deal.II/multigrid/mg_transfer_matrix_free.h

Lines 102 to 106 in 179a496

	void
	build(const DoFHandler<dim, dim> &dof_handler,
	const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
	&external_partitioners =
	std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>());

here (in that case, to avoid some vector copies between what kind of ghosts are needed for the level transfer and the matrix-vector product, to get another 5-10% of performance).

[peterrum]

@kronbichler I guess you know what the intention is for the new lamda function: the class MGTransferMatrixFree uses partitioners that are not compatible with the partioners used by MatrixFree. The solution for this problem appeared to me either to write adjust_ghost_range_if_necessary() or to overload copy_to_mg().

My question is the following: if I fill std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>> with the partitioners from MatrixFree do I get the same result as with the two approaches above? If yes, I don't see to add an forth way to accomplish the same thing! Furthermore if this is the case, I would like to make it also consistent in MGTransferGlobalCoarsening!

One minor question: why did you use std::vector and not MGLevelObject? I think I might have reviewed the patch that introduced this feature but I might have missed the intention back then!

[kronbichler]

Yes, the intention is to get the same result - the setup of the transfer class should check if the partitioner it needs is compatible with the one provided externally, which means that the MG transfer should only need ghosts which are already in the user-provided partitioner. In case it needs more indices, which happens sooner or later as we progress to coarser levels (in my experience, all levels but the finest one), it needs to keep a second set of local vectors and partitioners, which we then copy back and forth into. But since the problem sizes for all but the largest level are already ~8x smaller (in 3D and h-mg or p-mg with bisection), the copy for the small sizes is more tolerable as the largest one.

So we could use that approach in my opinion, but I think we still need to give it some thought as to which variant we prefer. I lean towards preferring your approach with lambdas, as the user does not need to collect things and can merely use the given level argument. Also, I have not checked yet if it obviates the adjust_ghost_range_if_necessary, i.e., if we set up the vectors with the correct partitioner. In case we do not, that would of course be the way to go.

[peterrum]

@kronbichler Thanks for the details! I will investigate this in the next days and will take a look if I can integrate the same optimizations in the global-coarsening approach.

Furthermore, I could image to move the lambda to the build() method and simply wrap the vector of partitioners inside of the lambda! For the time being, I let this PR as WIP.

[peterrum]

@kronbichler I have rebased this PR and have converted the vector of partitioners internally to a lambda function. What do you think. As far as I see it, the approach with adjust_ghost_range_if_necessary() is broken, since Multigrid will reset the partitioners in each cycle and I don't see a way to fix this. We might also want to consider to update step-59.

[kronbichler]

Do you mean that it calls adjust_ghost_range_if_necessary exactly once per level in each multigrid iteration (as the vector is reset upon copy_to_mg)? Either way, as discussed in the post below, I think we should improve the initialization of all vectors to avoid these kind of operations.

[peterrum]

Do you mean that it calls adjust_ghost_range_if_necessary exactly once per level in each multigrid iteration (as the vector is reset upon copy_to_mg)?

Yes. That is the case.

[kronbichler]

I have given this some more thought. I think experience has shown that we need a better way to initialize the vectors that we use inside the multigrid algorithm. However, let us discuss where the right place for this action, as I am not sure if the transfer classes are the right place:

• Eventually, things end up in the Multigrid class, so one would like to provide multigrid with the right set of vector for the solution, t and defect vectors, as these interact with the Matrix object and should not need to make the copies (aka adjust_ghost_values_if_necessary). Whether that is done by a lambda function as done here or a more traditional style by passing in a MGLevelObject<VectorType> & argument with appropriately initialized vectors, i.e., the right exemplar, might be up to debate. I traditionally opted for the latter, but a lambda as done here seems slightly favorable in terms of code compactness.
• The transfer classes sometimes need more ghost elements in the vectors than the mat-vec, but in case they do not, we do not want to copy vectors as that will double the timings on the finest level (mem-bandwidth limited code). Thus, we want to know inside the transfer whether the ghosts needed by the transfer (for LA::d::Vector::local_element()) are compatible with the ones in the rest of multigrid. We need this while setting things up, because we need to create the correct local numbering of ghosted entries, as we have learned in a few PRs this spring. That would indicate the transfer is appropriate anyway.
• If we can, we should be able to indicate when a vector will be overwritten later (e.g. solution in a typical V-cycle) or needs to be zeroed out explicitly (all but the last vector in local smoothing multigrid), so we can skip the zeroing if possible.

[bangerth]

Have you all given this issue some more thought?

[peterrum]

@kronbichler I have made the modification. Could you take a look at the last commit; it is not optimal, since I need somehow access to the partitioner, I need to create a temporal vector...

[kronbichler]

Suggested change

	New: The classes MGTransferMatrixFree::build() now also
	New: The class MGTransferMatrixFree::build() now also

[kronbichler]

Suggested change

	smoothers that are build around MatrixFree objects.
	smoothers that are built around MatrixFree objects.

[peterrum]

@kronbichler I have fixed both typos!

[peterrum]

/rebuild