46 #ifndef MUELU_MAXWELL_UTILS_DEF_HPP 47 #define MUELU_MAXWELL_UTILS_DEF_HPP 53 #include "Xpetra_Map.hpp" 54 #include "Xpetra_CrsMatrixUtils.hpp" 55 #include "Xpetra_MatrixUtils.hpp" 60 #include "MueLu_ThresholdAFilterFactory.hpp" 61 #include "MueLu_Utilities.hpp" 63 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 64 #include "MueLu_Utilities_kokkos.hpp" 68 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 70 #include <Thyra_VectorBase.hpp> 71 #include <Thyra_SolveSupportTypes.hpp> 73 #include <Stratimikos_DefaultLinearSolverBuilder.hpp> 75 #include "Teuchos_AbstractFactoryStd.hpp" 77 #ifdef HAVE_MUELU_IFPACK2 78 #include <Thyra_Ifpack2PreconditionerFactory.hpp> 85 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 RCP<Matrix> & D0_Matrix_,
97 Teuchos::ArrayRCP<bool> & BCrows_,
98 Teuchos::ArrayRCP<bool> & BCcols_,
99 Teuchos::ArrayRCP<bool> & BCdomain_,
100 bool & allEdgesBoundary_,
101 bool & allNodesBoundary_) {
106 int BCedgesLocal = 0;
107 int BCnodesLocal = 0;
108 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 119 auto BCrowsKokkos=BCrowsKokkos_;
120 Kokkos::parallel_reduce(BCrowsKokkos_.size(), KOKKOS_LAMBDA (
int i,
int & sum) {
125 auto BCdomainKokkos = BCdomainKokkos_;
126 Kokkos::parallel_reduce(BCdomainKokkos_.size(), KOKKOS_LAMBDA (
int i,
int & sum) {
127 if (BCdomainKokkos(i))
131 #endif // HAVE_MUELU_KOKKOS_REFACTOR 138 BCcols_.resize(D0_Matrix_->getColMap()->getNodeNumElements());
139 BCdomain_.resize(D0_Matrix_->getDomainMap()->getNodeNumElements());
142 for (
auto it = BCrows_.begin(); it != BCrows_.end(); ++it)
145 for (
auto it = BCdomain_.begin(); it != BCdomain_.end(); ++it)
150 MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCedgesLocal, BCedges_);
151 MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCnodesLocal, BCnodes_);
154 allEdgesBoundary_ = Teuchos::as<Xpetra::global_size_t>(BCedges_) >= D0_Matrix_->getRangeMap()->getGlobalNumElements();
155 allNodesBoundary_ = Teuchos::as<Xpetra::global_size_t>(BCnodes_) >= D0_Matrix_->getDomainMap()->getGlobalNumElements();
159 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
161 RCP<Matrix> & D0_Matrix_,
162 RCP<Matrix> & SM_Matrix_,
163 RCP<Matrix> & M1_Matrix_,
164 RCP<Matrix> & Ms_Matrix_) {
166 bool defaultFilter =
false;
171 if (parameterList_.get<
bool>(
"refmaxwell: filter D0",
true) && D0_Matrix_->getNodeMaxNumRowEntries()>2) {
175 fineLevel.
Set(
"A",D0_Matrix_);
176 fineLevel.
setlib(D0_Matrix_->getDomainMap()->lib());
179 fineLevel.
Request(
"A",ThreshFact.get());
180 ThreshFact->Build(fineLevel);
181 D0_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
184 defaultFilter =
true;
187 if (!M1_Matrix_.is_null() && parameterList_.get<
bool>(
"refmaxwell: filter M1", defaultFilter)) {
188 RCP<Vector> diag = VectorFactory::Build(M1_Matrix_->getRowMap());
190 M1_Matrix_->getLocalDiagCopy(*diag);
196 fineLevel.
Set(
"A",M1_Matrix_);
197 fineLevel.
setlib(M1_Matrix_->getDomainMap()->lib());
199 fineLevel.
Request(
"A",ThreshFact.get());
200 ThreshFact->Build(fineLevel);
201 M1_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
204 if (!Ms_Matrix_.is_null() && parameterList_.get<
bool>(
"refmaxwell: filter Ms", defaultFilter)) {
205 RCP<Vector> diag = VectorFactory::Build(Ms_Matrix_->getRowMap());
207 Ms_Matrix_->getLocalDiagCopy(*diag);
213 fineLevel.
Set(
"A",Ms_Matrix_);
214 fineLevel.
setlib(Ms_Matrix_->getDomainMap()->lib());
216 fineLevel.
Request(
"A",ThreshFact.get());
217 ThreshFact->Build(fineLevel);
218 Ms_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
221 if (!SM_Matrix_.is_null() && parameterList_.get<
bool>(
"refmaxwell: filter SM", defaultFilter)) {
222 RCP<Vector> diag = VectorFactory::Build(SM_Matrix_->getRowMap());
224 SM_Matrix_->getLocalDiagCopy(*diag);
230 fineLevel.
Set(
"A",SM_Matrix_);
231 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
233 fineLevel.
Request(
"A",ThreshFact.get());
234 ThreshFact->Build(fineLevel);
235 SM_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
242 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
245 RCP<const Import> xpImporter = A.getCrsGraph()->getImporter();
246 if (!xpImporter.is_null())
247 xpImporter->setDistributorParameters(matvecParams);
248 RCP<const Export> xpExporter = A.getCrsGraph()->getExporter();
249 if (!xpExporter.is_null())
250 xpExporter->setDistributorParameters(matvecParams);
256 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 258 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
259 struct StratimikosWrapperImpl {
260 static RCP<Thyra::PreconditionerBase<Scalar> > setupStratimikosPreconditioner(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A,
261 RCP<ParameterList> params) {
262 throw std::runtime_error(
"setupStratimikosPreconditioner: Requires Scalar=double");
266 template<
class LocalOrdinal,
class GlobalOrdinal,
class Node>
268 static RCP<Thyra::PreconditionerBase<double> > setupStratimikosPreconditioner(RCP<Xpetra::Matrix<double,LocalOrdinal,GlobalOrdinal,Node> > A,
269 RCP<ParameterList> params) {
273 RCP<const Thyra::LinearOpBase<Scalar> > thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(A)->getCrsMatrix());
275 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
276 typedef Thyra::PreconditionerFactoryBase<Scalar> Base;
277 typedef Thyra::MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> ImplMueLu;
278 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, ImplMueLu>(),
"MueLu");
279 #ifdef HAVE_MUELU_IFPACK2 281 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Impl;
282 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(),
"Ifpack2");
285 linearSolverBuilder.setParameterList(params);
290 auto precFactory = Thyra::createPreconditioningStrategy(linearSolverBuilder);
291 auto prec = precFactory->createPrec();
293 precFactory->initializePrec(Thyra::defaultLinearOpSource(thyraA), prec.get(), Thyra::SUPPORT_SOLVE_UNSPECIFIED);
300 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
301 RCP<Thyra::PreconditionerBase<Scalar> >
302 Maxwell_Utils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setupStratimikosPreconditioner(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A,
303 RCP<ParameterList> params) {
304 return StratimikosWrapperImpl<SC,LO,GO,NO>::setupStratimikosPreconditioner(A,params);
314 #define MUELU_MAXWELL_UTILS_SHORT 315 #endif //ifdef MUELU_MAXWELL_UTILS_DEF_HPP #define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
#define HAVE_MUELU_KOKKOS_REFACTOR
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
static void removeExplicitZeros(Teuchos::ParameterList ¶meterList, RCP< Matrix > &D0_Matrix, RCP< Matrix > &SM_Matrix, RCP< Matrix > &M1_Matrix, RCP< Matrix > &Ms_Matrix)
Remove explicit zeros.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
Namespace for MueLu classes and methods.
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
void setlib(Xpetra::UnderlyingLib lib2)
void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Kokkos::View< bool *, typename Node::device_type > &dirichletRows, Kokkos::View< bool *, typename Node::device_type > dirichletCols, Kokkos::View< bool *, typename Node::device_type > dirichletDomain)
Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
static void setMatvecParams(Matrix &A, RCP< ParameterList > matvecParams)
Sets matvec params on a matrix.
static void detectBoundaryConditionsSM(RCP< Matrix > &SM_Matrix, RCP< Matrix > &D0_Matrix, magnitudeType rowSumTol, int &BCedges, int &BCnodes, Teuchos::ArrayRCP< bool > &BCrows, Teuchos::ArrayRCP< bool > &BCcols, Teuchos::ArrayRCP< bool > &BCdomain, bool &allEdgesBoundary, bool &allNodesBoundary)
Detect Dirichlet boundary conditions.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomains)
void SetLevelID(int levelID)
Set level number.
void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, typename Node::device_type > &dirichletRows)
Factory for building a thresholded operator.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.