From 0930621c681c47144cb9893f089e59c7e94ae96a Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa.nebel@dr-blatt.de>
Date: Thu, 18 Jan 2024 11:55:43 +0100
Subject: [PATCH] Add test for building the AMG matrix hierarchy
This tests makes sure the matrices in the hierarchy get created in a deterministic way when the flag 'useFixedOrder' of the coarseningCriterion is set to true.
---
dune/istl/paamg/test/CMakeLists.txt | 5 ++
dune/istl/paamg/test/coarsentest.cc | 115 ++++++++++++++++++++++++++++
2 files changed, 120 insertions(+)
create mode 100644 dune/istl/paamg/test/coarsentest.cc
diff --git a/dune/istl/paamg/test/CMakeLists.txt b/dune/istl/paamg/test/CMakeLists.txt
index 43bc9bc3c..1b71b0fc7 100644
--- a/dune/istl/paamg/test/CMakeLists.txt
+++ b/dune/istl/paamg/test/CMakeLists.txt
@@ -92,6 +92,11 @@ dune_add_test(SOURCES hierarchytest.cc
TIMEOUT 600
CMAKE_GUARD MPI_FOUND)
+dune_add_test(SOURCES coarsentest.cc
+ MPI_RANKS 1 2 4
+ TIMEOUT 600
+ CMAKE_GUARD MPI_FOUND)
+
dune_add_test(NAME pamgtest
SOURCES parallelamgtest.cc
MPI_RANKS 1 2 4
diff --git a/dune/istl/paamg/test/coarsentest.cc b/dune/istl/paamg/test/coarsentest.cc
new file mode 100644
index 000000000..9f33fd81f
--- /dev/null
+++ b/dune/istl/paamg/test/coarsentest.cc
@@ -0,0 +1,115 @@
+// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
+// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <config.h>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/istl/paamg/matrixhierarchy.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+#include <dune/istl/schwarz.hh>
+#include "anisotropic.hh"
+
+constexpr int blockSize = 10;
+
+using Communication = Dune::OwnerOverlapCopyCommunication<int,int>;
+using ParallelIndexSet = Communication::ParallelIndexSet;
+using MatrixBlock = Dune::FieldMatrix<double,blockSize,blockSize>;
+using BCRSMat = Dune::BCRSMatrix<MatrixBlock>;
+using VectorBlock = Dune::FieldVector<double,blockSize>;
+using Vector = Dune::BlockVector<VectorBlock>;
+
+/**
+ * Check if the matrices A and B are the same.
+ *
+ * @param BCRSMat A
+ * @param BCRSMat B
+ * @return bool - True if A and B are the same, false if they are not.
+ */
+bool areSame(BCRSMat A, BCRSMat B)
+{
+ for (auto rowIt = A.begin(); rowIt != A.end(); ++rowIt)
+ for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+ for (size_t k = 0; k < blockSize; k++)
+ for (size_t l = 0; l < blockSize; l++) {
+ auto entryA = (*colIt)[k][l];
+ try {
+ auto entryB = B[rowIt.index()][colIt.index()][k][l];
+ } catch (Dune::ISTLError e) {
+ return false; // If the entry B[rowIt.index()][colIt.index()][k][l] does not exist, then the matrices are not the same.
+ }
+ if (std::abs((*colIt)[k][l]) > 1e-7 and std::abs(B[rowIt.index()][colIt.index()][k][l]) > 1e-7
+ and std::abs((*colIt)[k][l] - B[rowIt.index()][colIt.index()][k][l])/(*colIt)[k][l] > 1e-7) {
+ return false; // Check if the entries are approximately the same.
+ }
+ }
+ return true;
+}
+
+/**
+ * Builds two AMG matrix hierarchies and checks if they are equal.
+ * When setting up such matrix hierarchies on multiple processes, they get created in a non-deterministic way and the two matrix hierarchies might differ.
+ *
+ * @param int N Total number of unknowns
+ * @param int rank Rank of current process
+ * @param int noProcesses Number of all processes
+ * @return int - 0 if the hierarchies are the same, 1 if they are not.
+ */
+int testHierarchy(int N, int rank, int noProcesses)
+{
+ int n;
+ Communication pinfo(Dune::MPIHelper::getCommunicator());
+ ParallelIndexSet& indices = pinfo.indexSet();
+
+ using RemoteIndices = Dune::RemoteIndices<ParallelIndexSet>;
+ RemoteIndices& remoteIndices = pinfo.remoteIndices();
+
+ BCRSMat mat = setupAnisotropic2d<MatrixBlock>(N, indices, Dune::MPIHelper::getCommunication(), &n);
+
+ remoteIndices.rebuild<false>();
+
+ using OverlapFlags = Dune::NegateSet<Communication::OwnerSet>;
+ using Operator = Dune::OverlappingSchwarzOperator<BCRSMat,Vector,Vector,Communication>;
+ using Hierarchy = Dune::Amg::MatrixHierarchy<Operator,Communication>;
+
+ Operator op(mat, pinfo);
+ Hierarchy hierarchyA(
+ Dune::stackobject_to_shared_ptr(op),
+ Dune::stackobject_to_shared_ptr(pinfo));
+
+ Hierarchy hierarchyB(
+ Dune::stackobject_to_shared_ptr(op),
+ Dune::stackobject_to_shared_ptr(pinfo));
+
+ using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >;
+
+ Criterion criterionA(1000,10,1.2,1.6,Dune::Amg::AccumulationMode::successiveAccu,true);
+ Criterion criterionB(1000,10,1.2,1.6,Dune::Amg::AccumulationMode::successiveAccu,true);
+
+ hierarchyA.template build<OverlapFlags>(criterionA);
+ hierarchyB.template build<OverlapFlags>(criterionB);
+
+ MPI_Barrier(MPI_COMM_WORLD);
+ //Check if the second finest matrices are the same
+ auto finestMatrixA = hierarchyA.matrices().finest();
+ finestMatrixA++;
+ auto finestMatrixB = hierarchyB.matrices().finest();
+ finestMatrixB++;
+ if (!areSame(finestMatrixA->getmat(), finestMatrixB->getmat())) {
+ std::cerr << "During run with " << noProcesses << " processes: The matrix hierarchy on rank " << rank << " was created in a non-deterministic way." << std::endl;
+ return 1;
+ }
+ return 0;
+}
+
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ int N=10;
+
+ if(argc>1)
+ N = atoi(argv[1]);
+
+ return testHierarchy(N, mpiHelper.rank(), mpiHelper.size());
+}
--
GitLab