From a313dc7ef9f03a6f6185dc2aa6a0e854495bf44e Mon Sep 17 00:00:00 2001 From: rstumin Date: Tue, 23 Nov 2021 09:57:26 -0700 Subject: [PATCH] Added unit test for nullspace rotation computation --- .../MueLu_NullspaceFactory_def.hpp | 4 +- .../muelu/test/scaling/comp_rotations.xml | 59 +++++++ packages/muelu/test/unit_tests/CMakeLists.txt | 1 + .../test/unit_tests/NullspaceFactory.cpp | 146 ++++++++++++++++++ 4 files changed, 208 insertions(+), 2 deletions(-) create mode 100755 packages/muelu/test/scaling/comp_rotations.xml create mode 100644 packages/muelu/test/unit_tests/NullspaceFactory.cpp diff --git a/packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_NullspaceFactory_def.hpp b/packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_NullspaceFactory_def.hpp index f834fba7bcba..e235bf3d8624 100644 --- a/packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_NullspaceFactory_def.hpp +++ b/packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_NullspaceFactory_def.hpp @@ -190,7 +190,7 @@ namespace MueLu { nsValues[j*numPDEs + i] = 1.0; } } - if ( (int) nullspaceDim > numPDEs ) { + if (( (int) nullspaceDim > numPDEs ) && ((int) numPDEs > 1)) { /* xy rotation */ ArrayRCP nsValues = nullspace->getDataNonConst(numPDEs); int numBlocks = nsValues.size() / numPDEs; @@ -200,7 +200,7 @@ namespace MueLu { nsValues[j*numPDEs + 1] = (xvals[j]-cx); } } - if ( (int) nullspaceDim == numPDEs+3 ) { + if (( (int) nullspaceDim == numPDEs+3 ) && ((int) numPDEs > 2)) { /* xz rotation */ ArrayRCP nsValues = nullspace->getDataNonConst(numPDEs+1); int numBlocks = nsValues.size() / numPDEs; diff --git a/packages/muelu/test/scaling/comp_rotations.xml b/packages/muelu/test/scaling/comp_rotations.xml new file mode 100755 index 000000000000..a0be36f2ed29 --- /dev/null +++ b/packages/muelu/test/scaling/comp_rotations.xml @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/packages/muelu/test/unit_tests/CMakeLists.txt b/packages/muelu/test/unit_tests/CMakeLists.txt index 5c70b56165e9..95f071604c61 100644 --- a/packages/muelu/test/unit_tests/CMakeLists.txt +++ b/packages/muelu/test/unit_tests/CMakeLists.txt @@ -42,6 +42,7 @@ APPEND_SET(SOURCES MueLu_UnitTests.cpp MueLu_TestHelpers.cpp MultiVectorTransferFactory.cpp + NullspaceFactory.cpp ParameterList/FactoryFactory.cpp ParameterList/ParameterListInterpreter.cpp PermutedTransferFactory.cpp diff --git a/packages/muelu/test/unit_tests/NullspaceFactory.cpp b/packages/muelu/test/unit_tests/NullspaceFactory.cpp new file mode 100644 index 000000000000..fee460fb3227 --- /dev/null +++ b/packages/muelu/test/unit_tests/NullspaceFactory.cpp @@ -0,0 +1,146 @@ +// @HEADER +// +// *********************************************************************** +// +// MueLu: A package for multigrid based preconditioning +// Copyright 2012 Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact +// Jonathan Hu (jhu@sandia.gov) +// Andrey Prokopenko (aprokop@sandia.gov) +// Ray Tuminaro (rstumin@sandia.gov) +// +// *********************************************************************** +// +// @HEADER +#include "Teuchos_UnitTestHarness.hpp" +#include + +#include "MueLu_config.hpp" + +#include "MueLu_TestHelpers.hpp" +#include "MueLu_Version.hpp" + +#include +#include +#include "MueLu_Exceptions.hpp" +#include "MueLu_Utilities.hpp" +#include "MueLu_NullspaceFactory.hpp" +//#include "MueLu_FactoryManager.hpp" +#include "MueLu_CreateXpetraPreconditioner.hpp" + +namespace MueLuTests { + + + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(NullspaceFactory,3DElasticity, Scalar, LocalOrdinal, GlobalOrdinal, Node) + { +# include "MueLu_UseShortNames.hpp" + MUELU_TESTING_SET_OSTREAM; + MUELU_TESTING_LIMIT_SCOPE(Scalar,GlobalOrdinal,Node); + out << "version: " << MueLu::Version() << std::endl; + + using TST = Teuchos::ScalarTraits; + using magnitude_type = typename TST::magnitudeType; + using TMT = Teuchos::ScalarTraits; + using real_type = typename TST::coordinateType; + using RealValuedMultiVector = Xpetra::MultiVector; + + Level fineLevel, coarseLevel; + TestHelpers::TestFactory::createTwoLevelHierarchy(fineLevel, coarseLevel); + RCP > comm = Parameters::getDefaultComm(); + + Teuchos::ParameterList matrixList; + matrixList.set("nx", 2); + matrixList.set("ny", 3); + matrixList.set("nz", 4); + matrixList.set("mx", comm->getSize() ); + matrixList.set("my", 1); + matrixList.set("mz", 1); + matrixList.set("matrixType","Elasticity3D"); + RCP Op = TestHelpers::TestFactory::BuildMatrix(matrixList,TestHelpers::Parameters::getLib()); + Op->SetFixedBlockSize(3); + + fineLevel.Set("A",Op); + + Teuchos::ParameterList galeriList; + galeriList.set("nx", 2); + galeriList.set("ny", 3); + galeriList.set("nz", 4); + RCP fake = TestHelpers::TestFactory::Build1DPoisson(24); // need scalar PDE to produce proper coordinate array + RCP coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates("3D", fake->getRowMap(), galeriList); + fineLevel.Set("Coordinates", coordinates); + + RCP nspace = rcp(new NullspaceFactory()); + fineLevel.Request("Nullspace", &(*nspace)); + Teuchos::ParameterList nspaceList = *(nspace->GetValidParameterList()); + nspaceList.set("nullspace: calculate rotations", true); + nspace->SetParameterList(nspaceList); + fineLevel.Request(*nspace); + nspace->Build(fineLevel); + RCP nullSpace = fineLevel.Get >("Nullspace",nspace.get()); + + ArrayRCP xvals = coordinates->getDataNonConst(0); + ArrayRCP yvals = coordinates->getDataNonConst(1); + ArrayRCP zvals = coordinates->getDataNonConst(2); + RCP handCompNullSpace = MultiVectorFactory::Build(Op->getRowMap(), 6); + handCompNullSpace->putScalar(0.0); + ArrayRCP nsValues; + nsValues = handCompNullSpace->getDataNonConst(0); for (int j = 0; j < nsValues.size(); j +=3) nsValues[j] = 1.; + nsValues = handCompNullSpace->getDataNonConst(1); for (int j = 1; j < nsValues.size(); j +=3) nsValues[j] = 1.; + nsValues = handCompNullSpace->getDataNonConst(2); for (int j = 2; j < nsValues.size(); j +=3) nsValues[j] = 1.; + nsValues = handCompNullSpace->getDataNonConst(3); for (int j = 0; j < nsValues.size(); j +=3) nsValues[j] = -(yvals[j/3]-.5); + for (int j = 1; j < nsValues.size(); j +=3) nsValues[j] = (xvals[j/3]-.5); + nsValues = handCompNullSpace->getDataNonConst(4); for (int j = 1; j < nsValues.size(); j +=3) nsValues[j] = -(zvals[j/3]-.5); + for (int j = 2; j < nsValues.size(); j +=3) nsValues[j] = (yvals[j/3]-.5); + nsValues = handCompNullSpace->getDataNonConst(5); for (int j = 0; j < nsValues.size(); j +=3) nsValues[j] = -(zvals[j/3]-.5); + for (int j = 2; j < nsValues.size(); j +=3) nsValues[j] = (xvals[j/3]-.5); + + RCP diff = MultiVectorFactory::Build(Op->getRowMap(),6); + diff->putScalar(0.0); + //diff = fineNS + (-1.0)*(handCompNullSpace) + 0*diff + diff->update(1.0,*nullSpace,-1.0,*handCompNullSpace,0.0); + + Teuchos::Array norms(6); + diff->norm2(norms); + for (LocalOrdinal i=0; i<6; ++i) { + out << "||diff_" << i << "||_2 = " << norms[i] << std::endl; + TEST_EQUALITY(norms[i] < 100*TMT::eps(), true); + } + + } // NullSpace test + +#define MUELU_ETI_GROUP(Scalar, LO, GO, Node) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(NullspaceFactory,3DElasticity,Scalar,LO,GO,Node) +#include + +} // namespace MueLuTests