46 #ifndef THYRA_MUELU_TPETRA_Q2Q1PRECONDITIONER_FACTORY_DEF_HPP 47 #define THYRA_MUELU_TPETRA_Q2Q1PRECONDITIONER_FACTORY_DEF_HPP 49 #ifdef HAVE_MUELU_EXPERIMENTAL 53 #include <Thyra_DefaultPreconditioner.hpp> 54 #include <Thyra_TpetraLinearOp.hpp> 55 #include <Thyra_TpetraThyraWrappers.hpp> 57 #include <Teuchos_Ptr.hpp> 58 #include <Teuchos_TestForException.hpp> 59 #include <Teuchos_Assert.hpp> 60 #include <Teuchos_Time.hpp> 61 #include <Teuchos_FancyOStream.hpp> 62 #include <Teuchos_VerbosityLevel.hpp> 64 #include <Teko_Utilities.hpp> 75 #include "../../research/q2q1/MueLu_Q2Q1PFactory.hpp" 76 #include "../../research/q2q1/MueLu_Q2Q1uPFactory.hpp" 78 #include "MueLu_AmalgamationFactory.hpp" 80 #include "MueLu_BlockedDirectSolver.hpp" 81 #include "MueLu_BlockedPFactory.hpp" 82 #include "MueLu_BlockedRAPFactory.hpp" 83 #include "MueLu_BraessSarazinSmoother.hpp" 84 #include "MueLu_CoalesceDropFactory.hpp" 85 #include "MueLu_ConstraintFactory.hpp" 87 #include "MueLu_DirectSolver.hpp" 88 #include "MueLu_EminPFactory.hpp" 89 #include "MueLu_FactoryManager.hpp" 90 #include "MueLu_FilteredAFactory.hpp" 91 #include "MueLu_GenericRFactory.hpp" 93 #include "MueLu_PatternFactory.hpp" 94 #include "MueLu_SchurComplementFactory.hpp" 95 #include "MueLu_SmootherFactory.hpp" 96 #include "MueLu_SmootherPrototype.hpp" 97 #include "MueLu_SubBlockAFactory.hpp" 98 #include "MueLu_TpetraOperator.hpp" 99 #include "MueLu_TrilinosSmoother.hpp" 105 #define MUELU_GPD(name, type, defaultValue) \ 106 (paramList.isParameter(name) ? paramList.get<type>(name) : defaultValue) 110 using Teuchos::rcp_const_cast;
111 using Teuchos::rcp_dynamic_cast;
112 using Teuchos::ParameterList;
113 using Teuchos::ArrayView;
114 using Teuchos::ArrayRCP;
116 using Teuchos::Array;
119 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
124 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 typedef Thyra ::TpetraLinearOp<SC,LO,GO,NO> ThyraTpetraLinOp;
130 const RCP<const LinearOpBase<SC> > fwdOp = fwdOpSrc.getOp();
131 const RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = rcp_dynamic_cast<
const ThyraTpetraLinOp>(fwdOp);
132 const RCP<const TpetraLinOp> tpetraFwdOp = Teuchos::nonnull(thyraTpetraFwdOp) ? thyraTpetraFwdOp->getConstTpetraOperator() : Teuchos::null;
133 const RCP<const TpetraCrsMat> tpetraFwdCrsMat = rcp_dynamic_cast<
const TpetraCrsMat>(tpetraFwdOp);
135 return Teuchos::nonnull(tpetraFwdCrsMat);
138 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
139 RCP<PreconditionerBase<Scalar> >
141 return rcp(
new DefaultPreconditioner<SC>);
144 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
146 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> > &fwdOpSrc, PreconditionerBase<Scalar> *prec,
const ESupportSolveUse supportSolveUse)
const {
148 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
149 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
150 TEUCHOS_ASSERT(prec);
153 const RCP<const LinearOpBase<SC> > fwdOp = fwdOpSrc->getOp();
154 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
156 typedef Thyra::TpetraLinearOp<SC,LO,GO,NO> ThyraTpetraLinOp;
157 const RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = rcp_dynamic_cast<
const ThyraTpetraLinOp>(fwdOp);
158 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraTpetraFwdOp));
161 const RCP<const TpetraLinOp> tpetraFwdOp = thyraTpetraFwdOp->getConstTpetraOperator();
162 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdOp));
165 const RCP<const TpetraCrsMat> tpetraFwdCrsMat = rcp_dynamic_cast<
const TpetraCrsMat>(tpetraFwdOp);
166 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdCrsMat));
169 const Teuchos::Ptr<DefaultPreconditioner<SC> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<SC> *
>(prec));
170 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
173 const RCP<TpetraCrsMat> tpetraFwdCrsMatNonConst = rcp_const_cast<TpetraCrsMat>(tpetraFwdCrsMat);
178 ParameterList paramList = *paramList_;
181 RCP<MultiVector> coords, nullspace, velCoords, presCoords;
183 Teko::LinearOp thA11, thA12, thA21, thA11_9Pt;
184 if (paramList.isType<RCP<MultiVector> >(
"Coordinates")) { coords = paramList.get<RCP<MultiVector> >(
"Coordinates"); paramList.remove(
"Coordinates"); }
185 if (paramList.isType<RCP<MultiVector> >(
"Nullspace")) { nullspace = paramList.get<RCP<MultiVector> >(
"Nullspace"); paramList.remove(
"Nullspace"); }
186 if (paramList.isType<RCP<MultiVector> >(
"Velcoords")) { velCoords = paramList.get<RCP<MultiVector> >(
"Velcoords"); paramList.remove(
"Velcoords"); }
187 if (paramList.isType<RCP<MultiVector> >(
"Prescoords")) { presCoords = paramList.get<RCP<MultiVector> >(
"Prescoords"); paramList.remove(
"Prescoords"); }
188 if (paramList.isType<ArrayRCP<LO> > (
"p2vMap")) { p2vMap = paramList.get<ArrayRCP<LO> > (
"p2vMap"); paramList.remove(
"p2vMap"); }
189 if (paramList.isType<Teko::LinearOp> (
"A11")) { thA11 = paramList.get<Teko::LinearOp> (
"A11"); paramList.remove(
"A11"); }
190 if (paramList.isType<Teko::LinearOp> (
"A12")) { thA12 = paramList.get<Teko::LinearOp> (
"A12"); paramList.remove(
"A12"); }
191 if (paramList.isType<Teko::LinearOp> (
"A21")) { thA21 = paramList.get<Teko::LinearOp> (
"A21"); paramList.remove(
"A21"); }
192 if (paramList.isType<Teko::LinearOp> (
"A11_9Pt")) { thA11_9Pt = paramList.get<Teko::LinearOp> (
"A11_9Pt"); paramList.remove(
"A11_9Pt"); }
195 const RCP<MueLuOperator> mueluPrecOp = Q2Q1MkPrecond(paramList, velCoords, presCoords, p2vMap, thA11, thA12, thA21, thA11_9Pt);
197 const RCP<LinearOpBase<SC> > thyraPrecOp = Thyra::createLinearOp(RCP<TpetraLinOp>(mueluPrecOp));
198 defaultPrec->initializeUnspecified(thyraPrecOp);
201 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
203 uninitializePrec(PreconditionerBase<Scalar> *prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
205 TEUCHOS_ASSERT(prec);
208 const Teuchos::Ptr<DefaultPreconditioner<SC> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<SC> *
>(prec));
209 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
213 *fwdOp = Teuchos::null;
216 if (supportSolveUse) {
218 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
221 defaultPrec->uninitialize();
226 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
228 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
229 paramList_ = paramList;
232 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
239 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
242 RCP<ParameterList> savedParamList = paramList_;
243 paramList_ = Teuchos::null;
244 return savedParamList;
247 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
248 RCP<const ParameterList>
253 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
254 RCP<const ParameterList>
256 static RCP<const ParameterList> validPL;
258 if (validPL.is_null())
259 validPL = rcp(
new ParameterList());
264 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
265 RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
270 const ArrayRCP<LocalOrdinal>& p2vMap,
271 const Teko::LinearOp& thA11,
const Teko::LinearOp& thA12,
const Teko::LinearOp& thA21,
const Teko::LinearOp& thA11_9Pt)
const 293 const RCP<const Teuchos::Comm<int> > comm = velCoords->getMap()->getComm();
296 RCP<Thyra::LinearOpBase<SC> > ThNonConstA11 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11);
297 RCP<Thyra::LinearOpBase<SC> > ThNonConstA21 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA21);
298 RCP<Thyra::LinearOpBase<SC> > ThNonConstA12 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA12);
299 RCP<Thyra::LinearOpBase<SC> > ThNonConstA11_9Pt = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11_9Pt);
301 RCP<TP_Op> TpetA11 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11);
302 RCP<TP_Op> TpetA21 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA21);
303 RCP<TP_Op> TpetA12 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA12);
304 RCP<TP_Op> TpetA11_9Pt = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11_9Pt);
306 RCP<TP_Crs> TpetCrsA11 = rcp_dynamic_cast<TP_Crs>(TpetA11);
307 RCP<TP_Crs> TpetCrsA21 = rcp_dynamic_cast<TP_Crs>(TpetA21);
308 RCP<TP_Crs> TpetCrsA12 = rcp_dynamic_cast<TP_Crs>(TpetA12);
309 RCP<TP_Crs> TpetCrsA11_9Pt = rcp_dynamic_cast<TP_Crs>(TpetA11_9Pt);
321 RCP<const Map> domainMap2 = tmp_A_12->getDomainMap();
322 RCP<const Map> rangeMap2 = tmp_A_21->getRangeMap();
325 ArrayView<const GO> rangeElem2 = rangeMap2->getNodeElementList();
326 ArrayView<const GO> domainElem2 = domainMap2->getNodeElementList();
327 ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getNodeElementList();
328 ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getNodeElementList();
331 GO indexBase = domainMap2->getIndexBase();
333 Array<GO> newRowElem2(numRows2, 0);
335 newRowElem2[i] = numVel + rangeElem2[i];
337 RCP<const Map> newRangeMap2 = MapFactory::Build(lib, numRows2, newRowElem2, indexBase, comm);
340 Array<GO> newColElem2(numCols2, 0);
342 newColElem2[i] = numVel + domainElem2[i];
344 RCP<const Map> newDomainMap2 = MapFactory::Build(lib, numCols2, newColElem2, indexBase, comm);
346 RCP<Matrix> A_12 = MatrixFactory::Build(tmp_A_12->getRangeMap(), newDomainMap2, tmp_A_12->getNodeMaxNumRowEntries());
347 RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getNodeMaxNumRowEntries());
349 RCP<CrsMatrix> A_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_11) ->getCrsMatrix();
350 RCP<CrsMatrix> A_12_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_12) ->getCrsMatrix();
351 RCP<CrsMatrix> A_21_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_21) ->getCrsMatrix();
352 RCP<CrsMatrix> A_11_crs_9Pt = rcp_dynamic_cast<CrsMatrixWrap>(A_11_9Pt)->getCrsMatrix();
355 RCP<Matrix> A_22 = MatrixFactory::Build(newRangeMap2, newDomainMap2, 1);
356 RCP<CrsMatrix> A_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_22) ->getCrsMatrix();
359 Array<SC> smallVal(1, 1.0e-10);
364 ArrayView<const LO> inds;
365 ArrayView<const SC> vals;
366 for (
LO row = 0; row < as<LO>(numRows2); ++row) {
367 tmp_A_21->getLocalRowView(row, inds, vals);
369 size_t nnz = inds.size();
370 Array<GO> newInds(nnz, 0);
371 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
372 newInds[colID] = colElem1[inds[colID]];
374 A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
375 A_22_crs->insertGlobalValues(newRowElem2[row], Array<LO>(1, newRowElem2[row]), smallVal);
377 A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
378 A_22_crs->fillComplete(newDomainMap2, newRangeMap2);
380 RCP<Matrix> A_22 = Teuchos::null;
381 RCP<CrsMatrix> A_22_crs = Teuchos::null;
383 ArrayView<const LO> inds;
384 ArrayView<const SC> vals;
385 for (
LO row = 0; row < as<LO>(numRows2); ++row) {
386 tmp_A_21->getLocalRowView(row, inds, vals);
388 size_t nnz = inds.size();
389 Array<GO> newInds(nnz, 0);
390 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
391 newInds[colID] = colElem1[inds[colID]];
393 A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
395 A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
400 for (
LO row = 0; row < as<LO>(tmp_A_12->getRowMap()->getNodeNumElements()); ++row) {
401 tmp_A_12->getLocalRowView(row, inds, vals);
403 size_t nnz = inds.size();
404 Array<GO> newInds(nnz, 0);
405 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
406 newInds[colID] = newColElem2[inds[colID]];
408 A_12_crs->insertGlobalValues(rowElem1[row], newInds, vals);
410 A_12_crs->fillComplete(newDomainMap2, tmp_A_12->getRangeMap());
412 RCP<Matrix> A_12_abs = Absolute(*A_12);
413 RCP<Matrix> A_21_abs = Absolute(*A_21);
418 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
419 Teuchos::FancyOStream& out = *fancy;
420 out.setOutputToRootOnly(0);
424 SC dropTol = (paramList.get<
int>(
"useFilters") ? 0.06 : 0.00);
425 RCP<Matrix> filteredA = FilterMatrix(*A_11, *A_11, dropTol);
426 RCP<Matrix> filteredB = FilterMatrix(*BBt, *BBt_abs, dropTol);
428 RCP<Matrix> fA_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredA);
429 RCP<Matrix> fA_12_crs = Teuchos::null;
430 RCP<Matrix> fA_21_crs = Teuchos::null;
431 RCP<Matrix> fA_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredB);
434 std::vector<size_t> stridingInfo(1, 1);
435 int stridedBlockId = -1;
437 Array<GO> elementList(numVel+numPres);
438 Array<GO> velElem = A_12_crs->getRangeMap()->getNodeElementList();
439 Array<GO> presElem = A_21_crs->getRangeMap()->getNodeElementList();
443 RCP<const Map> fullMap = StridedMapFactory::Build(
Xpetra::UseTpetra, numVel+numPres, elementList(), indexBase, stridingInfo, comm);
445 std::vector<RCP<const Map> > partMaps(2);
446 partMaps[0] = StridedMapFactory::Build(
Xpetra::UseTpetra, numVel, velElem, indexBase, stridingInfo, comm);
447 partMaps[1] = StridedMapFactory::Build(
Xpetra::UseTpetra, numPres, presElem, indexBase, stridingInfo, comm, stridedBlockId, numVel);
450 RCP<const MapExtractor> mapExtractor = MapExtractorFactory::Build(fullMap, partMaps);
451 RCP<BlockedCrsMatrix> fA = rcp(
new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
452 fA->setMatrix(0, 0, fA_11_crs);
453 fA->setMatrix(0, 1, fA_12_crs);
454 fA->setMatrix(1, 0, fA_21_crs);
455 fA->setMatrix(1, 1, fA_22_crs);
462 SetDependencyTree(M, paramList);
464 RCP<Hierarchy> H = rcp(
new Hierarchy);
465 RCP<MueLu::Level> finestLevel = H->GetLevel(0);
466 finestLevel->Set(
"A", rcp_dynamic_cast<Matrix>(fA));
467 finestLevel->Set(
"p2vMap", p2vMap);
470 finestLevel->Set(
"AForPat", A_11_9Pt);
471 H->SetMaxCoarseSize(
MUELU_GPD(
"coarse: max size",
int, 1));
481 H->Keep(
"Ptent", M.
GetFactory(
"Ptent").get());
482 H->Setup(M, 0,
MUELU_GPD(
"max levels",
int, 3));
485 for (
int i = 1; i < H->GetNumLevels(); i++) {
486 RCP<Matrix> P = H->GetLevel(i)->template Get<RCP<Matrix> >(
"P");
487 RCP<BlockedCrsMatrix> Pcrs = rcp_dynamic_cast<BlockedCrsMatrix>(P);
488 RCP<Matrix> Pp = Pcrs->getMatrix(1,1);
489 RCP<Matrix> Pv = Pcrs->getMatrix(0,0);
499 std::string smootherType =
MUELU_GPD(
"smoother: type", std::string,
"vanka");
500 ParameterList smootherParams;
501 if (paramList.isSublist(
"smoother: params"))
502 smootherParams = paramList.sublist(
"smoother: params");
503 M.
SetFactory(
"Smoother", GetSmoother(smootherType, smootherParams,
false));
505 std::string coarseType =
MUELU_GPD(
"coarse: type", std::string,
"direct");
506 ParameterList coarseParams;
507 if (paramList.isSublist(
"coarse: params"))
508 coarseParams = paramList.sublist(
"coarse: params");
509 M.
SetFactory(
"CoarseSolver", GetSmoother(coarseType, coarseParams,
true));
511 #ifdef HAVE_MUELU_DEBUG 515 RCP<BlockedCrsMatrix> A = rcp(
new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
516 A->setMatrix(0, 0, A_11);
517 A->setMatrix(0, 1, A_12);
518 A->setMatrix(1, 0, A_21);
519 A->setMatrix(1, 1, A_22);
522 H->GetLevel(0)->Set(
"A", rcp_dynamic_cast<Matrix>(A));
524 H->Setup(M, 0, H->GetNumLevels());
529 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
530 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
540 RCP<GraphBase> filteredGraph;
546 level.
Set<RCP<Matrix> >(
"A", rcpFromRef(Pattern));
548 RCP<AmalgamationFactory> amalgFactory = rcp(
new AmalgamationFactory());
550 RCP<CoalesceDropFactory> dropFactory = rcp(
new CoalesceDropFactory());
551 ParameterList dropParams = *(dropFactory->GetValidParameterList());
552 dropParams.set(
"lightweight wrap",
true);
553 dropParams.set(
"aggregation: drop scheme",
"classical");
554 dropParams.set(
"aggregation: drop tol", dropTol);
556 dropFactory->SetParameterList(dropParams);
557 dropFactory->SetFactory(
"UnAmalgamationInfo", amalgFactory);
560 level.
Request(
"Graph", dropFactory.get());
561 dropFactory->Build(level);
563 level.
Get(
"Graph", filteredGraph, dropFactory.get());
566 RCP<Matrix> filteredA;
572 level.
Set(
"A", rcpFromRef(A));
573 level.
Set(
"Graph", filteredGraph);
574 level.
Set(
"Filtering",
true);
576 RCP<FilteredAFactory> filterFactory = rcp(
new FilteredAFactory());
577 ParameterList filterParams = *(filterFactory->GetValidParameterList());
580 filterParams.set(
"filtered matrix: reuse graph",
false);
581 filterParams.set(
"filtered matrix: use lumping",
false);
582 filterFactory->SetParameterList(filterParams);
585 level.
Request(
"A", filterFactory.get());
586 filterFactory->Build(level);
588 level.
Get(
"A", filteredA, filterFactory.get());
592 filteredA->resumeFill();
593 size_t numRows = filteredA->getRowMap()->getNodeNumElements();
594 for (
size_t i = 0; i < numRows; i++) {
595 ArrayView<const LO> inds;
596 ArrayView<const SC> vals;
597 filteredA->getLocalRowView(i, inds, vals);
599 size_t nnz = inds.size();
601 Array<SC> valsNew = vals;
604 SC diag = Teuchos::ScalarTraits<SC>::zero();
605 for (
size_t j = 0; j < inds.size(); j++) {
611 "No diagonal found");
615 valsNew[diagIndex] -= diag;
617 filteredA->replaceLocalValues(i, inds, valsNew);
619 filteredA->fillComplete();
624 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
638 RCP<FactoryManager> M11 = rcp(
new FactoryManager()), M22 = rcp(
new FactoryManager());
639 SetBlockDependencyTree(*M11, 0, 0,
"velocity", paramList);
640 SetBlockDependencyTree(*M22, 1, 1,
"pressure", paramList);
642 RCP<BlockedPFactory> PFact = rcp(
new BlockedPFactory());
643 ParameterList pParamList = *(PFact->GetValidParameterList());
644 pParamList.set(
"backwards",
true);
645 PFact->SetParameterList(pParamList);
646 PFact->AddFactoryManager(M11);
647 PFact->AddFactoryManager(M22);
650 RCP<GenericRFactory> RFact = rcp(
new GenericRFactory());
651 RFact->SetFactory(
"P", PFact);
654 RCP<MueLu::Factory > AcFact = rcp(
new BlockedRAPFactory());
655 AcFact->SetFactory(
"R", RFact);
656 AcFact->SetFactory(
"P", PFact);
662 RCP<MueLu::Factory> coarseFact = rcp(
new SmootherFactory(rcp(
new BlockedDirectSolver()), Teuchos::null));
667 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
675 typedef MueLu::Q2Q1PFactory <SC,LO,GO,NO> Q2Q1PFactory;
676 typedef MueLu::Q2Q1uPFactory <SC,LO,GO,NO> Q2Q1uPFactory;
679 RCP<SubBlockAFactory> AFact = rcp(
new SubBlockAFactory());
681 AFact->SetParameter(
"block row", Teuchos::ParameterEntry(row));
682 AFact->SetParameter(
"block col", Teuchos::ParameterEntry(col));
685 RCP<MueLu::Factory> Q2Q1Fact;
687 const bool isStructured =
false;
690 Q2Q1Fact = rcp(
new Q2Q1PFactory);
693 Q2Q1Fact = rcp(
new Q2Q1uPFactory);
694 ParameterList q2q1ParamList = *(Q2Q1Fact->GetValidParameterList());
695 q2q1ParamList.set(
"mode", mode);
696 if (paramList.isParameter(
"dump status"))
697 q2q1ParamList.set(
"dump status", paramList.get<
bool>(
"dump status"));
698 if (paramList.isParameter(
"phase2"))
699 q2q1ParamList.set(
"phase2", paramList.get<
bool>(
"phase2"));
700 Q2Q1Fact->SetParameterList(q2q1ParamList);
702 Q2Q1Fact->SetFactory(
"A", AFact);
705 RCP<PatternFactory> patternFact = rcp(
new PatternFactory);
706 ParameterList patternParams = *(patternFact->GetValidParameterList());
709 patternParams.set(
"emin: pattern order", 0);
710 patternFact->SetParameterList(patternParams);
711 patternFact->SetFactory(
"A", AFact);
712 patternFact->SetFactory(
"P", Q2Q1Fact);
715 RCP<ConstraintFactory> CFact = rcp(
new ConstraintFactory);
716 CFact->SetFactory(
"Ppattern", patternFact);
719 RCP<EminPFactory> EminPFact = rcp(
new EminPFactory());
720 ParameterList eminParams = *(EminPFact->GetValidParameterList());
721 if (paramList.isParameter(
"emin: num iterations"))
722 eminParams.set(
"emin: num iterations", paramList.get<
int>(
"emin: num iterations"));
723 if (mode ==
"pressure") {
724 eminParams.set(
"emin: iterative method",
"cg");
726 eminParams.set(
"emin: iterative method",
"gmres");
727 if (paramList.isParameter(
"emin: iterative method"))
728 eminParams.set(
"emin: iterative method", paramList.get<std::string>(
"emin: iterative method"));
730 EminPFact->SetParameterList(eminParams);
731 EminPFact->SetFactory(
"A", AFact);
732 EminPFact->SetFactory(
"Constraint", CFact);
733 EminPFact->SetFactory(
"P", Q2Q1Fact);
736 if (mode ==
"velocity" && (!paramList.isParameter(
"velocity: use transpose") || paramList.get<
bool>(
"velocity: use transpose") ==
false)) {
739 RCP<GenericRFactory> RFact = rcp(
new GenericRFactory());
740 RFact->SetFactory(
"P", EminPFact);
745 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
746 RCP<MueLu::FactoryBase>
748 GetSmoother(
const std::string& type,
const ParameterList& paramList,
bool coarseSolver)
const {
749 typedef Teuchos::ParameterEntry ParameterEntry;
760 RCP<SmootherPrototype> smootherPrototype;
761 if (type ==
"none") {
762 return Teuchos::null;
764 }
else if (type ==
"vanka") {
766 ParameterList schwarzList;
767 schwarzList.set(
"schwarz: overlap level", as<int>(0));
768 schwarzList.set(
"schwarz: zero starting solution",
false);
769 schwarzList.set(
"subdomain solver name",
"Block_Relaxation");
771 ParameterList& innerSolverList = schwarzList.sublist(
"subdomain solver parameters");
772 innerSolverList.set(
"partitioner: type",
"user");
773 innerSolverList.set(
"partitioner: overlap",
MUELU_GPD(
"partitioner: overlap",
int, 1));
774 innerSolverList.set(
"relaxation: type",
MUELU_GPD(
"relaxation: type", std::string,
"Gauss-Seidel"));
775 innerSolverList.set(
"relaxation: sweeps",
MUELU_GPD(
"relaxation: sweeps",
int, 1));
776 innerSolverList.set(
"relaxation: damping factor",
MUELU_GPD(
"relaxation: damping factor",
double, 0.5));
777 innerSolverList.set(
"relaxation: zero starting solution",
false);
780 std::string ifpackType =
"SCHWARZ";
782 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, schwarzList));
784 }
else if (type ==
"schwarz") {
786 std::string ifpackType =
"SCHWARZ";
788 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, paramList));
790 }
else if (type ==
"braess-sarazin") {
793 bool lumping =
MUELU_GPD(
"bs: lumping",
bool,
false);
795 RCP<SchurComplementFactory> schurFact = rcp(
new SchurComplementFactory());
796 schurFact->SetParameter(
"omega", ParameterEntry(omega));
797 schurFact->SetParameter(
"lumping", ParameterEntry(lumping));
801 RCP<SmootherPrototype> schurSmootherPrototype;
802 std::string schurSmootherType = (paramList.isParameter(
"schur smoother: type") ? paramList.get<std::string>(
"schur smoother: type") :
"RELAXATION");
803 if (schurSmootherType ==
"RELAXATION") {
804 ParameterList schurSmootherParams = paramList.sublist(
"schur smoother: params");
806 schurSmootherPrototype = rcp(
new TrilinosSmoother(schurSmootherType, schurSmootherParams));
808 schurSmootherPrototype = rcp(
new DirectSolver());
810 schurSmootherPrototype->SetFactory(
"A", schurFact);
812 RCP<SmootherFactory> schurSmootherFact = rcp(
new SmootherFactory(schurSmootherPrototype));
815 RCP<FactoryManager> braessManager = rcp(
new FactoryManager());
816 braessManager->SetFactory(
"A", schurFact);
817 braessManager->SetFactory(
"Smoother", schurSmootherFact);
818 braessManager->SetFactory(
"PreSmoother", schurSmootherFact);
819 braessManager->SetFactory(
"PostSmoother", schurSmootherFact);
820 braessManager->SetIgnoreUserData(
true);
822 smootherPrototype = rcp(
new BraessSarazinSmoother());
823 smootherPrototype->SetParameter(
"Sweeps", ParameterEntry(
MUELU_GPD(
"bs: sweeps",
int, 1)));
824 smootherPrototype->SetParameter(
"lumping", ParameterEntry(lumping));
825 smootherPrototype->SetParameter(
"Damping factor", ParameterEntry(omega));
826 smootherPrototype->SetParameter(
"q2q1 mode", ParameterEntry(
true));
827 rcp_dynamic_cast<BraessSarazinSmoother>(smootherPrototype)->AddFactoryManager(braessManager, 0);
829 }
else if (type ==
"ilu") {
830 std::string ifpackType =
"RILUK";
832 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, paramList));
834 }
else if (type ==
"direct") {
835 smootherPrototype = rcp(
new BlockedDirectSolver());
841 return coarseSolver ? rcp(
new SmootherFactory(smootherPrototype, Teuchos::null)) : rcp(
new SmootherFactory(smootherPrototype));
844 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
845 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
851 const CrsMatrixWrap& Awrap =
dynamic_cast<const CrsMatrixWrap&
>(A);
853 ArrayRCP<const size_t> iaA;
854 ArrayRCP<const LO> jaA;
855 ArrayRCP<const SC> valA;
856 Awrap.getCrsMatrix()->getAllValues(iaA, jaA, valA);
858 ArrayRCP<size_t> iaB (iaA .size());
859 ArrayRCP<LO> jaB (jaA .size());
860 ArrayRCP<SC> valB(valA.size());
861 for (
int i = 0; i < iaA .size(); i++) iaB [i] = iaA[i];
862 for (
int i = 0; i < jaA .size(); i++) jaB [i] = jaA[i];
863 for (
int i = 0; i < valA.size(); i++) valB[i] = Teuchos::ScalarTraits<SC>::magnitude(valA[i]);
866 RCP<CrsMatrix> Bcrs = rcp_dynamic_cast<CrsMatrixWrap>(B)->getCrsMatrix();
867 Bcrs->setAllValues(iaB, jaB, valB);
874 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
876 return "Thyra::MueLuTpetraQ2Q1PreconditionerFactory";
882 #endif // ifdef THYRA_MUELU_TPETRA_Q2Q1PRECONDITIONER_FACTORY_DEF_HPP Generic Smoother Factory for generating the smoothers of the MG hierarchy.
This class specifies the default factory that should generate some data on a Level if the data does n...
Teuchos::RCP< MueLu::TpetraOperator< SC, LO, GO, NO > > Q2Q1MkPrecond(const ParameterList ¶mList, const Teuchos::RCP< Tpetra::MultiVector< SC, LO, GO, NO > > &velCoords, const Teuchos::RCP< Tpetra::MultiVector< SC, LO, GO, NO > > &presCoords, const Teuchos::ArrayRCP< LO > &p2vMap, const Teko::LinearOp &thA11, const Teko::LinearOp &thA12, const Teko::LinearOp &thA21, const Teko::LinearOp &thA11_9Pt) const
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.
BraessSarazin smoother for 2x2 block matrices.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Teuchos::RCP< Xpetra::Matrix< SC, LO, GO, NO > > FilterMatrix(Xpetra::Matrix< SC, LO, GO, NO > &A, Xpetra::Matrix< SC, LO, GO, NO > &Pattern, SC dropTol) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
Factory for building blocked, segregated prolongation operators.
virtual const RCP< const Map > & getRowMap() const
bool isCompatible(const LinearOpSourceBase< SC > &fwdOp) const
Class that encapsulates external library smoothers.
Factory for building coarse matrices.
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
Base class for smoother prototypes.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
void SetBlockDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, LO row, LO col, const std::string &mode, const ParameterList ¶mList) const
Factory for building restriction operators using a prolongator factory.
virtual const RCP< const Map > & getColMap() const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Class that holds all level-specific information.
Factory for building a thresholded operator.
std::string description() const
AmalgamationFactory for subblocks of strided map based amalgamation data.
Factory for building the constraint operator.
Teuchos::RCP< PreconditionerBase< SC > > createPrec() const
direct solver for nxn blocked matrices
Teuchos::RCP< Xpetra::Matrix< SC, LO, GO, NO > > Absolute(const Xpetra::Matrix< SC, LO, GO, NO > &A) const
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Various adapters that will create a MueLu preconditioner that is a Tpetra::Operator.
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Set Factory.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Get factory associated with a particular data name.
MueLu representation of a graph.
Factory for building the Schur Complement for a 2x2 block matrix.
Factory for creating a graph base on a given matrix.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
RCP< MueLu::FactoryBase > GetSmoother(const std::string &type, const ParameterList ¶mList, bool coarseSolver) const
void SetLevelID(int levelID)
Set level number.
Factory for building nonzero patterns for energy minimization.
RCP< Xpetra::Matrix< SC, LO, GO, NO > > TpetraCrs_To_XpetraMatrix(const Teuchos::RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > &Atpetra)
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Factory for building Energy Minimization prolongators.
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
virtual Teuchos::RCP< const Map > getRangeMap() const=0
void SetDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, const ParameterList ¶mList) const
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
Exception throws to report errors in the internal logical of the program.
Factory for building filtered matrices using filtered graphs.
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< SC > > &fwdOp, PreconditionerBase< SC > *prec, const ESupportSolveUse supportSolveUse) const
#define MUELU_GPD(name, type, defaultValue)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual Teuchos::RCP< const Map > getDomainMap() const=0
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
MueLuTpetraQ2Q1PreconditionerFactory()
void uninitializePrec(PreconditionerBase< SC > *prec, Teuchos::RCP< const LinearOpSourceBase< SC > > *fwdOp, ESupportSolveUse *supportSolveUse) const
static const RCP< const NoFactory > getRCP()
Static Get() functions.