Skip to content

Commit

Permalink
MueLu: Adding filtering support for BrickAggregation
Browse files Browse the repository at this point in the history
  • Loading branch information
@csiefer2
csiefer2 committed Apr 14, 2022
1 parent b09a0a0 commit bc5684d
Show file tree
Hide file tree
Showing 3 changed files with 119 additions and 8 deletions.
13 changes: 9 additions & 4 deletions packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_decl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,18 +59,20 @@
#include "MueLu_SingleLevelFactoryBase.hpp"
#include "MueLu_BrickAggregationFactory_fwd.hpp"

#include "MueLu_GraphBase.hpp"
#include "MueLu_Graph_fwd.hpp"
#include "MueLu_LWGraph_fwd.hpp"
#include "MueLu_Level_fwd.hpp"
#include "MueLu_Aggregates_fwd.hpp"
#include "MueLu_Exceptions.hpp"
#include "MueLu_Utilities_fwd.hpp"

/*!
@class BrickAggregationFactory
@brief Aggregation method for generating "brick" aggregates.
@brief Aggregation method for generating "brick" aggregates. It also does "hotdogs" and "pancakes."
This factory can generate:
- In 2D, 2x2 and 3x3 aggregates
- In 3D, 2x2x2 or 3x3x3 aggregates
This factory can generate aggregates of size 1, 2 or 3 in each dimension, in any combination.
*/

namespace MueLu {
Expand Down Expand Up @@ -136,6 +138,9 @@ namespace MueLu {
void Setup(const RCP<const Teuchos::Comm<int> >& comm, const RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> >& coords, const RCP<const Map>& map) const;
RCP<container> Construct1DMap(const RCP<const Teuchos::Comm<int> >& comm, const ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& x) const;

void BuildGraph(Level& currentLevel, const RCP<Matrix>& A) const;


bool isDirichlet(LocalOrdinal LID) const;
bool isRoot (LocalOrdinal LID) const;
GlobalOrdinal getRoot (LocalOrdinal LID) const;
Expand Down
110 changes: 106 additions & 4 deletions packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,13 +61,15 @@
#include <Xpetra_MultiVector.hpp>
#include <Xpetra_MultiVectorFactory.hpp>


#include "MueLu_Aggregates.hpp"
#include "MueLu_Level.hpp"
#include "MueLu_MasterList.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_Utilities.hpp"
#include "MueLu_GraphBase.hpp"
#include "MueLu_Graph.hpp"
#include "MueLu_LWGraph.hpp"


namespace MueLu {

Expand All @@ -86,9 +88,6 @@ namespace MueLu {

validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory for matrix");
validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for coordinates");
validParamList->set< RCP<const FactoryBase> >("Graph" , Teuchos::null, "Generating factory of the graph");
validParamList->set< RCP<const FactoryBase> >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");

return validParamList;
}

Expand Down Expand Up @@ -192,6 +191,9 @@ namespace MueLu {
<< (nDim_ > 1 ? "x " + toString(by_) : "")
<< (nDim_ > 2 ? "x " + toString(bz_) : "") << std::endl;

// Build the graph
BuildGraph(currentLevel,A);

// Construct aggregates
RCP<Aggregates> aggregates = rcp(new Aggregates(colMap));
aggregates->setObjectLabel("Brick");
Expand Down Expand Up @@ -480,6 +482,106 @@ namespace MueLu {
}


template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildGraph(Level& currentLevel, const RCP<Matrix>& A) const {
// TODO: Currently only works w/ 1 DOF per node
double dirichletThreshold = 0.0;

if(bx_ > 1 && (nDim_ <= 1 || by_ > 1) && (nDim_ <=2 || bz_>1) ) {
FactoryMonitor m(*this, "Generating Graph (trivial)", currentLevel);
/*** Case 1: Use the matrix is the graph ***/
// Bricks are of non-trivial size in all active dimensions
RCP<GraphBase> graph = rcp(new Graph(A->getCrsGraph(), "graph of A"));
ArrayRCP<bool > boundaryNodes = Teuchos::arcp_const_cast<bool>(MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold));
graph->SetBoundaryNodeMap(boundaryNodes);

if (GetVerbLevel() & Statistics1) {
GO numLocalBoundaryNodes = 0;
GO numGlobalBoundaryNodes = 0;
for (LO i = 0; i < boundaryNodes.size(); ++i)
if (boundaryNodes[i])
numLocalBoundaryNodes++;
RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
}
Set(currentLevel, "DofsPerNode", 1);
Set(currentLevel, "Graph", graph);
Set(currentLevel, "Filtering",false);
}
else {
FactoryMonitor m(*this, "Generating Graph", currentLevel);
/*** Case 2: Dropping required ***/
// There is at least one active dimension in which we are not coarsening.
// Those connections need to be dropped
bool drop_x = (bx_ == 1);
bool drop_y = (nDim_> 1 && by_ == 1);
bool drop_z = (nDim_> 2 && bz_ == 1);

ArrayRCP<LO> rows (A->getLocalNumRows()+1);
ArrayRCP<LO> columns(A->getLocalNumEntries());

size_t N = A->getRowMap()->getLocalNumElements();

// FIXME: Do this on the host because indexing functions are host functions
auto G = A->getLocalMatrixHost().graph;
auto rowptr = G.row_map;
auto colind = G.entries;

int ct=0;
rows[0] = 0;
for(size_t row=0; row<N; row++) {
// NOTE: Assumes that the first part of the colmap is the rowmap
int ir,jr,kr;
LO row2 = A->getColMap()->getLocalElement(A->getRowMap()->getGlobalElement(row));
getIJK(row2,ir,jr,kr);

for(size_t cidx=rowptr[row]; cidx<rowptr[row+1]; cidx++) {
int ic,jc,kc;
LO col = colind[cidx];
getIJK(col,ic,jc,kc);

if( (row2 !=col) && ((drop_x && ir != ic) || (drop_y && jr != jc) || (drop_z && kr != kc) )) {
// Drop it
// printf("[%4d] DROP row = (%d,%d,%d) col = (%d,%d,%d)\n",(int)row,ir,jr,kr,ic,jc,kc);
}
else {
// Keep it
// printf("[%4d] KEEP row = (%d,%d,%d) col = (%d,%d,%d)\n",(int)row,ir,jr,kr,ic,jc,kc);
columns[ct] = col;
ct++;
}
}
rows[row+1] = ct;
}//end for

RCP<GraphBase> graph = rcp(new LWGraph(rows, columns, A->getRowMap(), A->getColMap(), "thresholded graph of A"));


ArrayRCP<bool > boundaryNodes = Teuchos::arcp_const_cast<bool>(MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold));
graph->SetBoundaryNodeMap(boundaryNodes);

if (GetVerbLevel() & Statistics1) {
GO numLocalBoundaryNodes = 0;
GO numGlobalBoundaryNodes = 0;
for (LO i = 0; i < boundaryNodes.size(); ++i)
if (boundaryNodes[i])
numLocalBoundaryNodes++;
RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
}
Set(currentLevel, "DofsPerNode", 1);
Set(currentLevel, "Graph", graph);
Set(currentLevel, "Filtering",true);
}//end else


}//end BuildGraph




} //namespace MueLu

#endif /* MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_ */
4 changes: 4 additions & 0 deletions packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1157,6 +1157,10 @@ namespace MueLu {
MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick z Dirichlet", bool, aggParams);
aggFactory->SetParameterList(aggParams);

// Unlike other factories, BrickAggregationFactory makes the Graph/DofsPerNode itself
manager.SetFactory("Graph", aggFactory);
manager.SetFactory("DofsPerNode", aggFactory);
manager.SetFactory("Filtering", aggFactory);
if (levelID > 1) {
// We check for levelID > 0, as in the interpreter aggFactory for
// levelID really corresponds to level 0. Managers are clunky, as they
Expand Down

0 comments on commit bc5684d

Please sign in to comment.