Skip to content

Commit

Permalink
Added unit test for nullspace rotation computation
Browse files Browse the repository at this point in the history
  • Loading branch information
@rstumin
rstumin committed Nov 23, 2021
1 parent 3fc6301 commit a313dc7
Show file tree
Hide file tree
Showing 4 changed files with 208 additions and 2 deletions.
4 changes: 2 additions & 2 deletions packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_NullspaceFactory_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Scalar> nsValues = nullspace->getDataNonConst(numPDEs);
int numBlocks = nsValues.size() / numPDEs;
Expand All @@ -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<Scalar> nsValues = nullspace->getDataNonConst(numPDEs+1);
int numBlocks = nsValues.size() / numPDEs;
Expand Down
59 changes: 59 additions & 0 deletions packages/muelu/test/scaling/comp_rotations.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
<ParameterList name="MueLu">

<!-- =========== GENERAL ================ -->
<Parameter name="problem: type" type="string" value="Elasticity-3D"/>
<Parameter name="verbosity" type="string" value="extreme"/>
<Parameter name="coarse: max size" type="int" value="10"/>
<Parameter name="multigrid algorithm" type="string" value="sa"/>

<!-- reduces setup cost for symmetric problems -->
<Parameter name="transpose: use implicit" type="bool" value="true"/>
<Parameter name="nullspace: calculate rotations" type="bool" value="true"/>

<!-- start of default values for general options (can be omitted) -->
<Parameter name="max levels" type="int" value="5"/>
<Parameter name="number of equations" type="int" value="3"/>
<Parameter name="sa: use filtered matrix" type="bool" value="true"/>
<!-- end of default values -->

<!-- =========== AGGREGATION =========== -->
<Parameter name="aggregation: type" type="string" value="uncoupled"/>
<Parameter name="aggregation: drop scheme" type="string" value="classical"/>
<Parameter name="aggregation: ordering" type="string" value="graph"/>
<!-- Uncomment the next line to enable dropping of weak connections, which can help AMG convergence
for anisotropic problems. The exact value is problem dependent. -->
<!-- <Parameter name="aggregation: drop tol" type="double" value="0.02"/> -->

<!-- =========== SMOOTHING =========== -->
<Parameter name="smoother: type" type="string" value="RELAXATION"/>
<ParameterList name="smoother: params">
<Parameter name="relaxation: type" type="string" value="Jacobi"/>
<Parameter name="relaxation: sweeps" type="int" value="2"/>
<Parameter name="relaxation: damping factor" type="double" value="0.7"/>
</ParameterList>

<!-- =========== REPARTITIONING =========== -->
<Parameter name="repartition: enable" type="bool" value="true"/>
<Parameter name="repartition: partitioner" type="string" value="zoltan2"/>
<Parameter name="repartition: start level" type="int" value="3"/>
<Parameter name="repartition: min rows per proc" type="int" value="800"/>
<Parameter name="repartition: max imbalance" type="double" value="1.1"/>
<Parameter name="repartition: remap parts" type="bool" value="false"/>
<!-- start of default values for repartitioning (can be omitted) -->
<Parameter name="repartition: remap parts" type="bool" value="true"/>
<Parameter name="repartition: rebalance P and R" type="bool" value="false"/>
<ParameterList name="repartition: params">
<Parameter name="algorithm" type="string" value="multijagged"/>
</ParameterList>

<ParameterList name="matrixmatrix: kernel params">
<Parameter name="compute global constants: temporaries" type="bool" value="false"/>
<Parameter name="compute global constants" type="bool" value="false"/>
</ParameterList>


<!-- end of default values -->

<Parameter name="use kokkos refactor" type="bool" value="false"/>

</ParameterList>
1 change: 1 addition & 0 deletions packages/muelu/test/unit_tests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
146 changes: 146 additions & 0 deletions packages/muelu/test/unit_tests/NullspaceFactory.cpp
View file
Original file line number Diff line number Diff line change
@@ -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 ([email protected])
// Andrey Prokopenko ([email protected])
// Ray Tuminaro ([email protected])
//
// ***********************************************************************
//
// @HEADER
#include "Teuchos_UnitTestHarness.hpp"
#include <Teuchos_ScalarTraits.hpp>

#include "MueLu_config.hpp"

#include "MueLu_TestHelpers.hpp"
#include "MueLu_Version.hpp"

#include <Xpetra_MultiVectorFactory.hpp>
#include <MueLu_VerbosityLevel.hpp>
#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<SC>;
using magnitude_type = typename TST::magnitudeType;
using TMT = Teuchos::ScalarTraits<magnitude_type>;
using real_type = typename TST::coordinateType;
using RealValuedMultiVector = Xpetra::MultiVector<real_type,LO,GO,NO>;

Level fineLevel, coarseLevel;
TestHelpers::TestFactory<Scalar, LO, GO, NO>::createTwoLevelHierarchy(fineLevel, coarseLevel);
RCP<const Teuchos::Comm<int> > 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<Matrix> Op = TestHelpers::TestFactory<Scalar, LO, GO, NO>::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<Matrix> fake = TestHelpers::TestFactory<Scalar, LO, GO, NO>::Build1DPoisson(24); // need scalar PDE to produce proper coordinate array
RCP<RealValuedMultiVector> coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<SC,LO,GO,Map,RealValuedMultiVector>("3D", fake->getRowMap(), galeriList);
fineLevel.Set("Coordinates", coordinates);

RCP<NullspaceFactory> 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<MultiVector> nullSpace = fineLevel.Get<RCP<MultiVector> >("Nullspace",nspace.get());

ArrayRCP<Scalar> xvals = coordinates->getDataNonConst(0);
ArrayRCP<Scalar> yvals = coordinates->getDataNonConst(1);
ArrayRCP<Scalar> zvals = coordinates->getDataNonConst(2);
RCP<MultiVector> handCompNullSpace = MultiVectorFactory::Build(Op->getRowMap(), 6);
handCompNullSpace->putScalar(0.0);
ArrayRCP<Scalar> 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<MultiVector> 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<typename TST::magnitudeType> 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 <MueLu_ETI_4arg.hpp>

} // namespace MueLuTests

0 comments on commit a313dc7

Please sign in to comment.