Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MueLu: adding example for structured aggregation and geometric interpolation #4170

Merged
merged 1 commit into from
Jan 11, 2019
Merged

MueLu: adding example for structured aggregation and geometric interpolation #4170

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 21 additions & 3 deletions packages/galeri/src-xpetra/Galeri_XpetraCartesian.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,20 +74,26 @@ namespace Galeri {
template <class LocalOrdinal, class GlobalOrdinal, class Map>
Teuchos::RCP<Map> Cartesian1D(const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
const GlobalOrdinal nx,
const GlobalOrdinal mx) {
const GlobalOrdinal mx,
Teuchos::ParameterList & list) {
if (nx <= 0 || mx <= 0 || (mx > nx))
throw Exception(__FILE__, __LINE__,
"Incorrect input parameter to Maps::Cartesian1D()",
"nx = " + toString(nx) +
", mx = " + toString(mx));

typedef GlobalOrdinal GO;
typedef LocalOrdinal LO;

int myPID = comm->getRank();

GO startx, endx;
Utils::getSubdomainData<GO>(nx, mx, myPID, startx, endx);

list.set("lnx", Teuchos::as<LO>(endx - startx));
list.set("lny", -1);
list.set("lnz", -1);

size_t numMyElements = endx - startx;
std::vector<GO> myGlobalElements(numMyElements);

Expand All @@ -104,7 +110,8 @@ namespace Galeri {
template <class LocalOrdinal, class GlobalOrdinal, class Map>
Teuchos::RCP<Map> Cartesian2D(const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
const GlobalOrdinal nx, const GlobalOrdinal ny,
const GlobalOrdinal mx, const GlobalOrdinal my) {
const GlobalOrdinal mx, const GlobalOrdinal my,
Teuchos::ParameterList & list) {
if (nx <= 0 || ny <= 0 || mx <= 0 || my <= 0 || (mx > nx) || (my > ny))
throw(Exception(__FILE__, __LINE__,
"Incorrect input parameter to Maps::Cartesian2D()",
Expand All @@ -114,13 +121,18 @@ namespace Galeri {
", my = " + toString(my)));

typedef GlobalOrdinal GO;
typedef LocalOrdinal LO;

int myPID = comm->getRank();

GO startx, starty, endx, endy;
Utils::getSubdomainData(nx, mx, myPID % mx, startx, endx);
Utils::getSubdomainData(ny, my, myPID / mx, starty, endy);

list.set("lnx", Teuchos::as<LO>(endx - startx));
list.set("lny", Teuchos::as<LO>(endy - starty));
list.set("lnz", -1);

size_t numMyElements = (endx - startx) * (endy - starty);
std::vector<GO> myGlobalElements(numMyElements);

Expand All @@ -138,7 +150,8 @@ namespace Galeri {
template <class LocalOrdinal, class GlobalOrdinal, class Map>
Teuchos::RCP<Map> Cartesian3D(const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
const GlobalOrdinal nx, const GlobalOrdinal ny, const GlobalOrdinal nz,
const GlobalOrdinal mx, const GlobalOrdinal my, const GlobalOrdinal mz) {
const GlobalOrdinal mx, const GlobalOrdinal my, const GlobalOrdinal mz,
Teuchos::ParameterList & list) {
if (nx <= 0 || ny <= 0 || nz <= 0 ||
mx <= 0 || my <= 0 || mz <= 0 ||
(mx > nx) || (my > ny) || (mz > nz))
Expand All @@ -152,6 +165,7 @@ namespace Galeri {
", mz = " + toString(mz));

typedef GlobalOrdinal GO;
typedef LocalOrdinal LO;

GO mxy = mx * my;

Expand All @@ -162,6 +176,10 @@ namespace Galeri {
Utils::getSubdomainData(ny, my, (myPID % mxy) / mx, starty, endy);
Utils::getSubdomainData(nz, mz, myPID / mxy , startz, endz);

list.set("lnx", Teuchos::as<LO>(endx - startx));
list.set("lny", Teuchos::as<LO>(endy - starty));
list.set("lnz", Teuchos::as<LO>(endz - startz));

size_t numMyElements = (endx - startx) * (endy - starty) * (endz - startz);
std::vector<GO> myGlobalElements(numMyElements);

Expand Down
6 changes: 3 additions & 3 deletions packages/galeri/src-xpetra/Galeri_XpetraMaps.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,7 @@ namespace Galeri {
list.set("mx", mx);
}

return Maps::Cartesian1D<LocalOrdinal, GlobalOrdinal, Map>(comm, nx, mx);
return Maps::Cartesian1D<LocalOrdinal, GlobalOrdinal, Map>(comm, nx, mx, list);

} else if (mapType == "Cartesian2D") {
if (nx == -1 || ny == -1) {
Expand Down Expand Up @@ -253,7 +253,7 @@ namespace Galeri {
list.set("my", my);
}

return Maps::Cartesian2D<LocalOrdinal, GlobalOrdinal, Map>(comm, nx, ny, mx, my);
return Maps::Cartesian2D<LocalOrdinal, GlobalOrdinal, Map>(comm, nx, ny, mx, my, list);

} else if (mapType == "Cartesian3D") {
if (nx == -1 || ny == -1 || nz == -1) {
Expand Down Expand Up @@ -318,7 +318,7 @@ namespace Galeri {
+ ", mz = " + toString(mz));
}

return Maps::Cartesian3D<LocalOrdinal, GlobalOrdinal, Map>(comm, nx, ny, nz, mx, my, mz);
return Maps::Cartesian3D<LocalOrdinal, GlobalOrdinal, Map>(comm, nx, ny, nz, mx, my, mz, list);

} else {
throw Exception(__FILE__, __LINE__,
Expand Down
8 changes: 7 additions & 1 deletion packages/muelu/src/Graph/StructuredAggregation/MueLu_AggregationStructuredAlgorithm_def.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,7 @@ namespace MueLu {
void AggregationStructuredAlgorithm<LocalOrdinal, GlobalOrdinal, Node>::
BuildGraph(const GraphBase& graph, RCP<IndexManager>& geoData, RCP<CrsGraph>& myGraph,
RCP<const Map>& coarseCoordinatesFineMap, RCP<const Map>& coarseCoordinatesMap) const {
Monitor m(*this, "BuildAggregates");
Monitor m(*this, "BuildGraphP");

RCP<Teuchos::FancyOStream> out;
if(const char* dbg = std::getenv("MUELU_STRUCTUREDALGORITHM_DEBUG")) {
Expand Down Expand Up @@ -242,14 +242,20 @@ namespace MueLu {
graph.GetDomainMap()->getNode());
}

*out << "Call constructor of CrsGraph" << std::endl;
myGraph = CrsGraphFactory::Build(graph.GetDomainMap(),
colMap,
nnzOnRow,
Xpetra::DynamicProfile);

*out << "Fill CrsGraph" << std::endl;
for(LO nodeIdx = 0; nodeIdx < geoData->getNumLocalFineNodes(); ++nodeIdx) {
myGraph->insertLocalIndices(nodeIdx, colIndex(rowPtr[nodeIdx], nnzOnRow[nodeIdx]) );
}

*out << "Call fillComplete on CrsGraph" << std::endl;
myGraph->fillComplete(domainMap, graph.GetDomainMap());
*out << "Prolongator CrsGraph computed" << std::endl;

} // BuildAggregates()

Expand Down
22 changes: 18 additions & 4 deletions packages/muelu/src/Graph/StructuredAggregation/MueLu_StructuredAggregationFactory_def.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,10 +100,12 @@ namespace MueLu {

validParamList->set<RCP<const FactoryBase> >("Graph", Teuchos::null,
"Graph of the matrix after amalgamation but without dropping.");
validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
"Number of spatial dimension provided by CoordinatesTransferFactory.");
validParamList->set<RCP<const FactoryBase> >("gNodesPerDim", Teuchos::null,
"Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
"Global number of nodes per spatial dimension provided by CoordinatesTransferFactory.");
validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
"Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
"Local number of nodes per spatial dimension provided by CoordinatesTransferFactory.");

return validParamList;
} // GetValidParameterList
Expand Down Expand Up @@ -133,6 +135,13 @@ namespace MueLu {

// Request the local number of nodes per dimensions
if(currentLevel.GetLevelID() == 0) {
if(currentLevel.IsAvailable("numDimensions", NoFactory::get())) {
currentLevel.DeclareInput("numDimensions", NoFactory::get(), this);
} else {
TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("numDimensions", NoFactory::get()),
Exceptions::RuntimeError,
"numDimensions was not provided by the user on level0!");
}
if(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
currentLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
} else {
Expand All @@ -141,6 +150,7 @@ namespace MueLu {
"lNodesPerDim was not provided by the user on level0!");
}
} else {
Input(currentLevel, "numDimensions");
Input(currentLevel, "lNodesPerDim");
}
} // DeclareInput
Expand Down Expand Up @@ -172,23 +182,25 @@ namespace MueLu {

// Since we want to operate on nodes and not dof, we need to modify the rowMap in order to
// obtain a nodeMap.
const int numDimensions = pL.get<int>("aggregation: number of spatial dimensions");
const int interpolationOrder = pL.get<int>("aggregation: coarsening order");
std::string meshLayout = pL.get<std::string>("aggregation: mesh layout");
std::string coupling = pL.get<std::string>("aggregation: coupling");
const bool coupled = (coupling == "coupled" ? true : false);
std::string outputType = pL.get<std::string>("aggregation: output type");
const bool outputAggregates = (outputType == "Aggregates" ? true : false);
int numDimensions;
Array<GO> gFineNodesPerDir(3);
Array<LO> lFineNodesPerDir(3);
if(currentLevel.GetLevelID() == 0) {
// On level 0, data is provided by applications and has no associated factory.
numDimensions = currentLevel.Get<int>("numDimensions", NoFactory::get());
if(coupled) {
gFineNodesPerDir = currentLevel.Get<Array<GO> >("gNodesPerDim", NoFactory::get());
}
lFineNodesPerDir = currentLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
} else {
// On level > 0, data is provided directly by generating factories.
numDimensions = Get<int>(currentLevel, "numDimensions");
if(coupled) {
gFineNodesPerDir = Get<Array<GO> >(currentLevel, "gNodesPerDim");
}
Expand Down Expand Up @@ -274,6 +286,8 @@ namespace MueLu {


*out << "The index manager has now been built" << std::endl;
*out << "graph num nodes: " << fineMap->getNodeNumElements()
<< ", structured aggregation num nodes: " << geoData->getNumLocalFineNodes() << std::endl;
TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
!= static_cast<size_t>(geoData->getNumLocalFineNodes()),
Exceptions::RuntimeError,
Expand Down Expand Up @@ -330,7 +344,7 @@ namespace MueLu {
Set(currentLevel, "coarseCoordinatesFineMap", coarseCoordinatesFineMap);
Set(currentLevel, "coarseCoordinatesMap", coarseCoordinatesMap);
Set(currentLevel, "interpolationOrder", interpolationOrder);
Set(currentLevel, "numDimensions", numDimensions);
// Set(currentLevel, "coarseNumDimensions", numDimensions);

} // Build
} //namespace MueLu
Expand Down
5 changes: 5 additions & 0 deletions packages/muelu/src/Misc/MueLu_CoordinatesTransferFactory_def.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ namespace MueLu {
validParamList->set<RCP<const FactoryBase> >("coarseCoordinates", Teuchos::null, "Factory for coarse coordinates generation");
validParamList->set<RCP<const FactoryBase> >("gCoarseNodesPerDim", Teuchos::null, "Factory providing the global number of nodes per spatial dimensions of the mesh");
validParamList->set<RCP<const FactoryBase> >("lCoarseNodesPerDim", Teuchos::null, "Factory providing the local number of nodes per spatial dimensions of the mesh");
validParamList->set<RCP<const FactoryBase> >("numDimensions" , Teuchos::null, "Factory providing the number of spatial dimensions of the mesh");
validParamList->set<int> ("write start", -1, "first level at which coordinates should be written to file");
validParamList->set<int> ("write end", -1, "last level at which coordinates should be written to file");
validParamList->set<RCP<const FactoryBase> >("aggregationRegionTypeCoarse", Teuchos::null, "Factory indicating what aggregation type is to be used on the coarse level of the region");
Expand All @@ -93,6 +94,7 @@ namespace MueLu {
Input(fineLevel, "gCoarseNodesPerDim");
}
Input(fineLevel, "lCoarseNodesPerDim");
Input(fineLevel, "numDimensions");
} else if(pL.get<bool>("Geometric") == true) {
Input(coarseLevel, "coarseCoordinates");
Input(coarseLevel, "gCoarseNodesPerDim");
Expand Down Expand Up @@ -121,6 +123,7 @@ namespace MueLu {

GetOStream(Runtime0) << "Transferring coordinates" << std::endl;

int numDimensions;
RCP<xdMV> coarseCoords;
RCP<xdMV> fineCoords;
Array<GO> gCoarseNodesPerDir;
Expand All @@ -142,6 +145,8 @@ namespace MueLu {
}
lCoarseNodesPerDir = Get<Array<LO> >(fineLevel, "lCoarseNodesPerDim");
Set<Array<LO> >(coarseLevel, "lNodesPerDim", lCoarseNodesPerDir);
numDimensions = Get<int>(fineLevel, "numDimensions");
Set<int>(coarseLevel, "numDimensions", numDimensions);
} else if(pL.get<bool>("Geometric") == true) {
coarseCoords = Get<RCP<xdMV> >(coarseLevel, "coarseCoordinates");
gCoarseNodesPerDir = Get<Array<GO> >(coarseLevel, "gCoarseNodesPerDim");
Expand Down
Loading