47 #ifndef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP 48 #define THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP 52 #ifdef HAVE_MUELU_STRATIMIKOS 55 #include "Thyra_DefaultPreconditioner.hpp" 56 #include "Thyra_BlockedLinearOpBase.hpp" 57 #include "Thyra_DiagonalLinearOpBase.hpp" 59 #ifdef HAVE_MUELU_TPETRA 60 #include "Thyra_TpetraLinearOp.hpp" 61 #include "Thyra_TpetraThyraWrappers.hpp" 63 #ifdef HAVE_MUELU_EPETRA 64 #include "Thyra_EpetraLinearOp.hpp" 65 #include "Thyra_EpetraThyraWrappers.hpp" 68 #include "Teuchos_Ptr.hpp" 70 #include "Teuchos_Assert.hpp" 73 #include <Xpetra_CrsMatrixWrap.hpp> 74 #include <Xpetra_CrsMatrix.hpp> 75 #include <Xpetra_Matrix.hpp> 76 #include <Xpetra_ThyraUtils.hpp> 78 #include <MueLu_Hierarchy.hpp> 80 #include <MueLu_HierarchyUtils.hpp> 81 #include <MueLu_Utilities.hpp> 82 #include <MueLu_ParameterListInterpreter.hpp> 83 #include <MueLu_MLParameterListInterpreter.hpp> 86 #include <MueLu_RefMaxwell.hpp> 87 #ifdef HAVE_MUELU_TPETRA 88 #include <MueLu_TpetraOperator.hpp> 90 #ifdef HAVE_MUELU_EPETRA 92 #include <Xpetra_EpetraOperator.hpp> 95 #include "Thyra_PreconditionerFactoryBase.hpp" 97 #include "Kokkos_DefaultNode.hpp" 110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node = KokkosClassic::DefaultNode::DefaultNodeType>
125 bool isCompatible(
const LinearOpSourceBase<Scalar>& fwdOp)
const;
130 PreconditionerBase<Scalar>* prec,
131 const ESupportSolveUse supportSolveUse
136 ESupportSolveUse* supportSolveUse
172 #ifdef HAVE_MUELU_EPETRA 203 #ifdef HAVE_MUELU_TPETRA 204 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
207 #ifdef HAVE_MUELU_EPETRA 208 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp))
return true;
211 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp))
return true;
223 PreconditionerBase<Scalar>* prec,
224 const ESupportSolveUse supportSolveUse
226 using Teuchos::rcp_dynamic_cast;
229 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
230 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
231 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
232 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
233 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
234 typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
235 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
236 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
237 typedef Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
238 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
239 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
240 #ifdef HAVE_MUELU_TPETRA 242 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 243 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 244 typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
247 #if defined(HAVE_MUELU_EPETRA) 248 typedef Thyra::EpetraLinearOp ThyEpLinOp;
249 typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> XpEpCrsMat;
269 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
270 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
271 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
279 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
291 RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
295 RCP<XpMat> A00 =
rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
312 RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
316 A =
rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
326 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
334 const bool startingOver = (thyra_precOp.
is_null() || !paramList.
isParameter(
"reuse: type") || paramList.
get<std::string>(
"reuse: type") ==
"none");
336 if (startingOver ==
true) {
343 #ifdef HAVE_MUELU_TPETRA 345 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 346 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 352 paramList.
remove(
"Nullspace");
353 RCP<XpMultVec> nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_nullspace);
370 RCP<XpCrsMat> crsM1NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM1);
373 RCP<XpMat> M1 =
rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM1NonConst));
394 RCP<XpCrsMat> crsD0NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsD0);
397 RCP<XpMat> D0 =
rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsD0NonConst));
409 paramList.
remove(
"M0inv");
410 RCP<XpCrsMat> xM0inv = rcp_dynamic_cast<XpCrsMat>(tM0inv,
true);
414 paramList.
remove(
"M0inv");
417 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
418 RCP<XpMat> M0inv = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(Xpetra::toXpetra(tDiag));
422 paramList.
remove(
"M0inv");
426 RCP<XpCrsMat> crsM0invNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM0inv);
429 RCP<XpMat> M0inv =
rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM0invNonConst));
439 "Thyra::MueLuRefMaxwellPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
443 #ifdef HAVE_MUELU_EPETRA 448 paramList.
remove(
"Nullspace");
451 RCP<XpMultVec> nullspace = rcp_dynamic_cast<XpMultVec>(xpEpNullspaceMult,
true);
460 RCP<XpCrsMat> xCrsM1 = rcp_dynamic_cast<XpCrsMat>(xeM1,
true);
462 RCP<XpMat> xM1 = rcp_dynamic_cast<XpMat>(xwM1);
471 RCP<XpCrsMat> crsM1NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM1);
474 RCP<XpMat> M1 =
rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM1NonConst));
488 RCP<XpCrsMat> xCrsD0 = rcp_dynamic_cast<XpCrsMat>(xeD0,
true);
490 RCP<XpMat> xD0 = rcp_dynamic_cast<XpMat>(xwD0);
499 RCP<XpCrsMat> crsD0NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsD0);
502 RCP<XpMat> D0 =
rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsD0NonConst));
514 paramList.
remove(
"M0inv");
516 RCP<XpCrsMat> xCrsM0inv = rcp_dynamic_cast<XpCrsMat>(xeM0inv,
true);
518 RCP<XpMat> xM0inv = rcp_dynamic_cast<XpMat>(xwM0inv);
523 paramList.
remove(
"M0inv");
533 RCP<XpMat> M0inv = Xpetra::MatrixFactory<double,int,int,Node>::Build(xpDiag);
537 paramList.
remove(
"M0inv");
541 RCP<XpCrsMat> crsM0invNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM0inv);
544 RCP<XpMat> M0inv =
rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM0invNonConst));
555 paramList.
set<
bool>(
"refmaxwell: use as preconditioner",
true);
562 #if defined(HAVE_MUELU_TPETRA) 564 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 565 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 566 RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
570 "Thyra::MueLuRefMaxwellPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
574 #if defined(HAVE_MUELU_EPETRA)// && defined(HAVE_MUELU_SERIAL) 576 RCP<ThyEpLinOp> epetr_precOp = rcp_dynamic_cast<ThyEpLinOp>(thyra_precOp);
589 RCP<XpOp> xpOp = Teuchos::rcp_dynamic_cast<XpOp>(preconditioner);
590 thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
594 defaultPrec->initializeUnspecified(thyraPrecOp);
600 ESupportSolveUse* supportSolveUse
610 *fwdOp = Teuchos::null;
613 if (supportSolveUse) {
615 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
618 defaultPrec->uninitialize();
635 return savedParamList;
656 std::string
description()
const {
return "Thyra::MueLuRefMaxwellPreconditionerFactory"; }
666 #endif // HAVE_MUELU_EPETRA 670 #endif // #ifdef HAVE_MUELU_STRATIMIKOS 672 #endif // THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
T & get(const std::string &name, T def_value)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool nonnull(const std::shared_ptr< T > &p)
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
MueLuRefMaxwellPreconditionerFactory()
bool is_null(const std::shared_ptr< T > &p)
Concrete preconditioner factory subclass for Thyra based on MueLu.Add support for MueLu preconditione...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::string description() const
Teuchos::RCP< Teuchos::ParameterList > paramList_
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
static RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList ¶mList)
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOpSrc) const
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form...
std::string description() const
bool remove(std::string const &name, bool throwIfNotExists=true)
MueLuRefMaxwellPreconditionerFactory()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< Teuchos::ParameterList > paramList_
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOp, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
bool isType(const std::string &name) const
Teuchos::RCP< PreconditionerBase< Scalar > > createPrec() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
Teuchos::RCP< PreconditionerBase< Scalar > > createPrec() const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOp) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
bool isParameter(const std::string &name) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Exception throws to report errors in the internal logical of the program.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)