Compare MGTransferMatrixFree and MGTransferMF for compute_no_normal_flux_constraints_on_level()

[peterrum]

This PR shows that MGTransferMatrixFree and MGTransferMF have different behaviors in the case that DoFs are constrained by other DoFs (as in the case of no-normal-flux constrains). The difference is that MGTransferMF zeroes the constrained values during prolongation while MGTransferMatrixFree does not. The wouldn't be a problem (if one just ignores these values in the vector). However, restriction is the transposed operation of prolongation; in the case that some values are in the fine vector, this leads to different restricted residuals.

In #16167, I could observe that, in the case that constrains are not zeroed, the results are different if the system matrix has zeroes or ones on the diagonal entries related to constrained DoFs.

I noticed this discrepancy while compare times of MGTransferMatrixFree and MGTransferMF. I noticed that

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

Lines 2237 to 2286 in 087064a

	// ------------------------------- weights -------------------------------
	if (transfer.fine_element_is_continuous)
	{
	// compute weights globally
	LinearAlgebra::distributed::Vector<Number> weight_vector;
	compute_weights(constraints_fine, transfer, weight_vector);
	// ... and store them cell-wise a DG format
	transfer.weights.resize(n_dof_indices_fine.back());
	Number *weights_0 = transfer.weights.data() + n_dof_indices_fine[0];
	Number *weights_1 = transfer.weights.data() + n_dof_indices_fine[1];
	unsigned int *dof_indices_fine_0 =
	transfer.constraint_info_fine.dof_indices.data() +
	n_dof_indices_fine[0];
	unsigned int *dof_indices_fine_1 =
	transfer.constraint_info_fine.dof_indices.data() +
	n_dof_indices_fine[1];
	process_cells(
	[&](const auto &, const auto &) {
	for (unsigned int i = 0;
	i < transfer.schemes[0].n_dofs_per_cell_fine;
	++i)
	weights_0[i] =
	weight_vector.local_element(dof_indices_fine_0[i]);
	dof_indices_fine_0 += transfer.schemes[0].n_dofs_per_cell_fine;
	weights_0 += transfer.schemes[0].n_dofs_per_cell_fine;
	},
	[&](const auto &, const auto &, const auto c) {
	for (unsigned int i = 0;
	i < transfer.schemes[1].n_dofs_per_cell_coarse;
	++i)
	weights_1[cell_local_children_indices[c][i]] =
	weight_vector.local_element(
	dof_indices_fine_1[cell_local_children_indices[c][i]]);
	// move pointers (only once at the end)
	if (c + 1 == GeometryInfo<dim>::max_children_per_cell)
	{
	dof_indices_fine_1 +=
	transfer.schemes[1].n_dofs_per_cell_fine;
	weights_1 += transfer.schemes[1].n_dofs_per_cell_fine;
	}
	});
	if (is_feq)
	compress_weights(transfer);
	}

add additional costs due to additional work.

Q: Since we want to merge both classes. We need to decide which is the correct behavior and make the (incompatible) change now.

[kronbichler]

Q: Since we want to merge both classes. We need to decide which is the correct behavior and make the (incompatible) change now.

I agree, we should make this change now and unify the code.

However, restriction is the transposed operation of prolongation; in the case that some values are in the fine vector, this leads to different restricted residuals.

Please help me to better understand the problem: We get different values in case the restrict_and_add function call if the input vector contains non-zero entries in constrained DoFs. We did not see this problem before because we usually have the right-hand side and residual zero for constrained DoFs. Does that sound right?

Regarding the desired behavior, I would suggest to go with the more robust behavior in terms of solvers. Does it even make sense to compute something for the constrained DoFs, or would any such value likely be problematic because it might get added with other contributions that then need to be addressed by smoothers? My question is motivated by the fact that the transfer matrices are rectangular, so I considered rank-deficient matrices and similar artifacts the common case. In other words, my philosophy was that we should keep all non-zero contributions out of the mg cycle, which would mean that we should ignore the vector entries in constrained DoFs. I could be wrong, and one might in fact be able to have useful information on constrained DoFs, but then we would have to resolve the constraints properly, wouldn't we?

[peterrum]

We get different values in case the restrict_and_add function call if the input vector contains non-zero entries in constrained DoFs.

Exactly.

We did not see this problem before because we usually have the right-hand side and residual zero for constrained DoFs. Does that sound right?

Yes. At least all tests in multigrid-global-coarsening did this by using MatrixFree (and as a result implicitly zeroing the DoFs). However, this is an assumption we cannot make; see for example https://github.com/exadg/exadg/blob/56272c15cceb9ecc1300a176c7c5f11cb78604ca/include/exadg/operators/operator_base.cpp#L270-L274 (which was the motivation for #16167).

Regarding the desired behavior, I would suggest to go with the more robust behavior in terms of solvers. Does it even make sense to compute something for the constrained DoFs, or would any such value likely be problematic because it might get added with other contributions that then need to be addressed by smoothers? My question is motivated by the fact that the transfer matrices are rectangular, so I considered rank-deficient matrices and similar artifacts the common case. In other words, my philosophy was that we should keep all non-zero contributions out of the mg cycle, which would mean that we should ignore the vector entries in constrained DoFs. I could be wrong, and one might in fact be able to have useful information on constrained DoFs, but then we would have to resolve the constraints properly, wouldn't we?

Exactly, this is the reasoning for zeroing constrained values during weighting. However, this is not done in the case of MGTransferMatrixFree.

It would be a bit odd that MGTransferMatrixFree only takes care of the constraints on the coarse level but not on the fine level with the result that the results differ if the level operator has ones or zeroes on the diagonal.

[kronbichler]

I agree, we should not consider constraints on either side.

A longer-term desire would be that we should have a way to get rid of all constrained DoFs while solving linear systems with a matrix-free solver (to avoid any kind of ambiguities, because those cases are purely artificial). We would need some preparations for making that work, in particular a way for MatrixFree to initialize vectors compared to a DoFHandler that carries more DoFs, and mechanisms to exchange between the vectors. It will not solve the problem here, but it would probably avoid a series of issues in application codes if we had a way to do it. See #15224.

[peterrum]

As a note, MGConstrainedDoFs::make_no_normal_flux_constraints() only supports the case that the normal are aligned with x-, y-, or z-axis, with the result that the components are not constrained with each other and one can treat the constrain as DBC on a single component. Maybe this is done due to known limitations?

To apply more complex constraints (e.g., compute_no_normal_flux_constraints_on_level(); see for example https://github.com/geodynamics/aspect/blob/4fc9e5ee44e5cc66b4583d82b5c446eb2b7e486a/source/simulator/stokes_matrix_free.cc#L2519C1-L2529 in ASPECT), one needs to use MGConstrainedDoFs::add_user_constraints(). But this seems to be very fragile and unsafe, since arbitrary constrains can be applied without any internal checks (e.g., inhomogeneous DBCs). These features are hardly tested in deal.II.

[tjhei]

Maybe this is done due to known limitations?

This is done because MGConstrainedDoFs only supports 0 constraints (efficiently) and custom constraints (inefficiently). We want to use the former if possible.
We could have made make_no_normal_flux_constraints put things into the user constraints if they are not axis aligned. Does it make sense to make that change?

These features are hardly tested in deal.II.

Do you mean VectorTools::compute_no_normal_flux_constraints_on_level() or using it in Multigrid?

[peterrum]

@kronbichler With the last PR, I apply the constraints consistently on the fine level (but also on the coarse level during restriction via condense()).

[kronbichler]

I guess this one can be closed now?