MueLu: Cut Drop Boundary Nodes Not Behaving As Intended

[NyquilDreams]

Bug Report

@jhux2 @cgcgcg @csiefer2

Description

While testing the Cut Drop Algorithm in MueLu_CoalesceDropFactory_def.hpp using the unit tests in CoalesceDropFactory.cpp, I noticed that the unit tests did not test for a matrix with boundary nodes. I then noticed that the boundaryNodes View in the algorithm was erroneously being accessed using colID (a loop index) instead of col (the local column index).

Trilinos/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp

Lines 555 to 567 in bf922e7

	for (LO colID = 0; colID < (LO)nnz; colID++) {
	LO col = indices[colID];
	if (row == col) {
	drop_vec.emplace_back(zero, one, colID, false);
	continue;
	}
	// Don't aggregate boundaries
	if (boundaryNodes[colID]) continue;
	typename STS::magnitudeType aiiajj = STS::magnitude(threshold * threshold * ghostedDiagVals[col] * ghostedDiagVals[row]); // eps^2*|a_ii|*|a_jj|
	typename STS::magnitudeType aij = STS::magnitude(vals[colID] * vals[colID]); // |a_ij|^2
	drop_vec.emplace_back(aij, aiiajj, colID, false);
	}

However, an incorrect assumption is also being made. Boundary nodes is obtained using DetectDirichletRows_kokkos_host which returns a View with a size equal to the number of local rows. This does not account for overlap.

Trilinos/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp

Line 458 in bf922e7

auto boundaryNodes = MueLu::Utilities<SC, LO, GO, NO>::DetectDirichletRows_kokkos_host(*A, dirichletThreshold);

This seems to have been noticed further down in the code. The boundaryNodes View used by the Cut Drop algorithm should make similar changes to the following excerpt. The unit test should also be modified accordingly.

Trilinos/packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp

Lines 1868 to 1876 in bf922e7

	// ghosting of boundary node map
	auto boundaryNodes = inputGraph->GetBoundaryNodeMap();
	auto boundaryNodesVector = Xpetra::VectorFactory<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>::Build(inputGraph->GetDomainMap());
	for (size_t i = 0; i < inputGraph->GetNodeNumVertices(); i++)
	boundaryNodesVector->getDataNonConst(0)[i] = boundaryNodes[i];
	// Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write("boundary",*boundaryNodesVector);
	auto boundaryColumnVector = Xpetra::VectorFactory<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>::Build(inputGraph->GetImportMap());
	boundaryColumnVector->doImport(*boundaryNodesVector, *importer, Xpetra::INSERT);
	auto boundaryColumn = boundaryColumnVector->getData(0);

[github-actions]

Automatic mention of the @trilinos/muelu team