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>
66#include <Xpetra_BlockedCrsMatrix.hpp>
67#include <Xpetra_CrsMatrixWrap.hpp>
68#include <Xpetra_IO.hpp>
69#include <Xpetra_MapExtractorFactory.hpp>
70#include <Xpetra_Matrix.hpp>
71#include <Xpetra_MatrixMatrix.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;
127 typedef Tpetra::Operator <SC,LO,GO,NO> TpetraLinOp;
128 typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TpetraCrsMat;
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));
160 typedef Tpetra::Operator<SC,LO,GO,NO> TpetraLinOp;
161 const RCP<const TpetraLinOp> tpetraFwdOp = thyraTpetraFwdOp->getConstTpetraOperator();
162 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdOp));
164 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraCrsMat;
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_;
180 typedef Tpetra::MultiVector<SC,LO,GO,NO> MultiVector;
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> >
268 const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& velCoords,
269 const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& presCoords,
270 const ArrayRCP<LocalOrdinal>& p2vMap,
271 const Teko::LinearOp& thA11,
const Teko::LinearOp& thA12,
const Teko::LinearOp& thA21,
const Teko::LinearOp& thA11_9Pt)
const
275 typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TP_Crs;
276 typedef Tpetra::Operator <SC,LO,GO,NO> TP_Op;
278 typedef Xpetra::BlockedCrsMatrix <SC,LO,GO,NO> BlockedCrsMatrix;
279 typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
280 typedef Xpetra::CrsMatrixWrap <SC,LO,GO,NO> CrsMatrixWrap;
281 typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
282 typedef Xpetra::MapExtractor <SC,LO,GO,NO> MapExtractor;
283 typedef Xpetra::Map <LO,GO,NO> Map;
284 typedef Xpetra::MapFactory <LO,GO,NO> MapFactory;
285 typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
286 typedef Xpetra::MatrixFactory <SC,LO,GO,NO> MatrixFactory;
287 typedef Xpetra::StridedMapFactory <LO,GO,NO> StridedMapFactory;
289 typedef MueLu::Hierarchy <SC,LO,GO,NO> Hierarchy;
291 const RCP<const Teuchos::Comm<int> > comm = velCoords->getMap()->getComm();
294 RCP<Thyra::LinearOpBase<SC> > ThNonConstA11 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11);
295 RCP<Thyra::LinearOpBase<SC> > ThNonConstA21 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA21);
296 RCP<Thyra::LinearOpBase<SC> > ThNonConstA12 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA12);
297 RCP<Thyra::LinearOpBase<SC> > ThNonConstA11_9Pt = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11_9Pt);
299 RCP<TP_Op> TpetA11 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11);
300 RCP<TP_Op> TpetA21 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA21);
301 RCP<TP_Op> TpetA12 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA12);
302 RCP<TP_Op> TpetA11_9Pt = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11_9Pt);
304 RCP<TP_Crs> TpetCrsA11 = rcp_dynamic_cast<TP_Crs>(TpetA11);
305 RCP<TP_Crs> TpetCrsA21 = rcp_dynamic_cast<TP_Crs>(TpetA21);
306 RCP<TP_Crs> TpetCrsA12 = rcp_dynamic_cast<TP_Crs>(TpetA12);
307 RCP<TP_Crs> TpetCrsA11_9Pt = rcp_dynamic_cast<TP_Crs>(TpetA11_9Pt);
314 Xpetra::global_size_t numVel = A_11->getRowMap()->getLocalNumElements();
315 Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getLocalNumElements();
319 RCP<const Map> domainMap2 = tmp_A_12->getDomainMap();
320 RCP<const Map> rangeMap2 = tmp_A_21->getRangeMap();
321 Xpetra::global_size_t numRows2 = rangeMap2->getLocalNumElements();
322 Xpetra::global_size_t numCols2 = domainMap2->getLocalNumElements();
323 ArrayView<const GO> rangeElem2 = rangeMap2->getLocalElementList();
324 ArrayView<const GO> domainElem2 = domainMap2->getLocalElementList();
325 ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getLocalElementList();
326 ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getLocalElementList();
328 Xpetra::UnderlyingLib lib = domainMap2->lib();
329 GO indexBase = domainMap2->getIndexBase();
331 Array<GO> newRowElem2(numRows2, 0);
332 for (Xpetra::global_size_t i = 0; i < numRows2; i++)
333 newRowElem2[i] = numVel + rangeElem2[i];
335 RCP<const Map> newRangeMap2 = MapFactory::Build(lib, numRows2, newRowElem2, indexBase, comm);
338 Array<GO> newColElem2(numCols2, 0);
339 for (Xpetra::global_size_t i = 0; i < numCols2; i++)
340 newColElem2[i] = numVel + domainElem2[i];
342 RCP<const Map> newDomainMap2 = MapFactory::Build(lib, numCols2, newColElem2, indexBase, comm);
344 RCP<Matrix> A_12 = MatrixFactory::Build(tmp_A_12->getRangeMap(), newDomainMap2, tmp_A_12->getLocalMaxNumRowEntries());
345 RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getLocalMaxNumRowEntries());
347 RCP<CrsMatrix> A_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_11) ->getCrsMatrix();
348 RCP<CrsMatrix> A_12_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_12) ->getCrsMatrix();
349 RCP<CrsMatrix> A_21_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_21) ->getCrsMatrix();
350 RCP<CrsMatrix> A_11_crs_9Pt = rcp_dynamic_cast<CrsMatrixWrap>(A_11_9Pt)->getCrsMatrix();
353 RCP<Matrix> A_22 = MatrixFactory::Build(newRangeMap2, newDomainMap2, 1);
354 RCP<CrsMatrix> A_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_22) ->getCrsMatrix();
357 Array<SC> smallVal(1, 1.0e-10);
362 ArrayView<const LO> inds;
363 ArrayView<const SC> vals;
364 for (
LO row = 0; row < as<LO>(numRows2); ++row) {
365 tmp_A_21->getLocalRowView(row, inds, vals);
367 size_t nnz = inds.size();
368 Array<GO> newInds(nnz, 0);
369 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
370 newInds[colID] = colElem1[inds[colID]];
372 A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
373 A_22_crs->insertGlobalValues(newRowElem2[row], Array<LO>(1, newRowElem2[row]), smallVal);
375 A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
376 A_22_crs->fillComplete(newDomainMap2, newRangeMap2);
378 RCP<Matrix> A_22 = Teuchos::null;
379 RCP<CrsMatrix> A_22_crs = Teuchos::null;
381 ArrayView<const LO> inds;
382 ArrayView<const SC> vals;
383 for (
LO row = 0; row < as<LO>(numRows2); ++row) {
384 tmp_A_21->getLocalRowView(row, inds, vals);
386 size_t nnz = inds.size();
387 Array<GO> newInds(nnz, 0);
388 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
389 newInds[colID] = colElem1[inds[colID]];
391 A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
393 A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
398 for (
LO row = 0; row < as<LO>(tmp_A_12->getRowMap()->getLocalNumElements()); ++row) {
399 tmp_A_12->getLocalRowView(row, inds, vals);
401 size_t nnz = inds.size();
402 Array<GO> newInds(nnz, 0);
403 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
404 newInds[colID] = newColElem2[inds[colID]];
406 A_12_crs->insertGlobalValues(rowElem1[row], newInds, vals);
408 A_12_crs->fillComplete(newDomainMap2, tmp_A_12->getRangeMap());
410 RCP<Matrix> A_12_abs = Absolute(*A_12);
411 RCP<Matrix> A_21_abs = Absolute(*A_21);
416 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
417 Teuchos::FancyOStream& out = *fancy;
418 out.setOutputToRootOnly(0);
419 RCP<Matrix> BBt = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21,
false, *A_12,
false, out);
420 RCP<Matrix> BBt_abs = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21_abs,
false, *A_12_abs,
false, out);
422 SC dropTol = (paramList.get<
int>(
"useFilters") ? paramList.get<
double>(
"tau_1") : 0.00);
423 RCP<Matrix> filteredA = FilterMatrix(*A_11, *A_11, dropTol);
424 RCP<Matrix> filteredB = FilterMatrix(*BBt, *BBt_abs, dropTol);
426 RCP<Matrix> fA_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredA);
427 RCP<Matrix> fA_12_crs = Teuchos::null;
428 RCP<Matrix> fA_21_crs = Teuchos::null;
429 RCP<Matrix> fA_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredB);
432 std::vector<size_t> stridingInfo(1, 1);
433 int stridedBlockId = -1;
435 Array<GO> elementList(numVel+numPres);
436 Array<GO> velElem = A_12_crs->getRangeMap()->getLocalElementList();
437 Array<GO> presElem = A_21_crs->getRangeMap()->getLocalElementList();
439 for (Xpetra::global_size_t i = 0 ; i < numVel; i++) elementList[i] = velElem[i];
440 for (Xpetra::global_size_t i = numVel; i < numVel+numPres; i++) elementList[i] = presElem[i-numVel];
441 RCP<const Map> fullMap = StridedMapFactory::Build(Xpetra::UseTpetra, numVel+numPres, elementList(), indexBase, stridingInfo, comm);
443 std::vector<RCP<const Map> > partMaps(2);
444 partMaps[0] = StridedMapFactory::Build(Xpetra::UseTpetra, numVel, velElem, indexBase, stridingInfo, comm);
445 partMaps[1] = StridedMapFactory::Build(Xpetra::UseTpetra, numPres, presElem, indexBase, stridingInfo, comm, stridedBlockId, numVel);
448 RCP<const MapExtractor> mapExtractor = MapExtractorFactory::Build(fullMap, partMaps);
449 RCP<BlockedCrsMatrix> fA = rcp(
new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
450 fA->setMatrix(0, 0, fA_11_crs);
451 fA->setMatrix(0, 1, fA_12_crs);
452 fA->setMatrix(1, 0, fA_21_crs);
453 fA->setMatrix(1, 1, fA_22_crs);
460 SetDependencyTree(M, paramList);
462 RCP<Hierarchy> H = rcp(
new Hierarchy);
463 RCP<MueLu::Level> finestLevel = H->GetLevel(0);
464 finestLevel->Set(
"A", rcp_dynamic_cast<Matrix>(fA));
465 finestLevel->Set(
"p2vMap", p2vMap);
466 finestLevel->Set(
"CoordinatesVelocity", Xpetra::toXpetra(velCoords));
467 finestLevel->Set(
"CoordinatesPressure", Xpetra::toXpetra(presCoords));
468 finestLevel->Set(
"AForPat", A_11_9Pt);
469 H->SetMaxCoarseSize(
MUELU_GPD(
"coarse: max size",
int, 1));
479 H->Keep(
"Ptent", M.
GetFactory(
"Ptent").get());
480 H->Setup(M, 0,
MUELU_GPD(
"max levels",
int, 3));
483 for (
int i = 1; i < H->GetNumLevels(); i++) {
484 RCP<Matrix> P = H->GetLevel(i)->template Get<RCP<Matrix> >(
"P");
485 RCP<BlockedCrsMatrix> Pcrs = rcp_dynamic_cast<BlockedCrsMatrix>(P);
486 RCP<Matrix> Pp = Pcrs->getMatrix(1,1);
487 RCP<Matrix> Pv = Pcrs->getMatrix(0,0);
489 Xpetra::IO<SC,LO,GO,NO>::Write(
"Pp_l" +
MueLu::toString(i) +
".mm", *Pp);
490 Xpetra::IO<SC,LO,GO,NO>::Write(
"Pv_l" +
MueLu::toString(i) +
".mm", *Pv);
497 std::string smootherType =
MUELU_GPD(
"smoother: type", std::string,
"vanka");
498 ParameterList smootherParams;
499 if (paramList.isSublist(
"smoother: params"))
500 smootherParams = paramList.sublist(
"smoother: params");
501 M.
SetFactory(
"Smoother", GetSmoother(smootherType, smootherParams,
false));
503 std::string coarseType =
MUELU_GPD(
"coarse: type", std::string,
"direct");
504 ParameterList coarseParams;
505 if (paramList.isSublist(
"coarse: params"))
506 coarseParams = paramList.sublist(
"coarse: params");
507 M.
SetFactory(
"CoarseSolver", GetSmoother(coarseType, coarseParams,
true));
509#ifdef HAVE_MUELU_DEBUG
513 RCP<BlockedCrsMatrix> A = rcp(
new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
514 A->setMatrix(0, 0, A_11);
515 A->setMatrix(0, 1, A_12);
516 A->setMatrix(1, 0, A_21);
517 A->setMatrix(1, 1, A_22);
520 H->GetLevel(0)->Set(
"A", rcp_dynamic_cast<Matrix>(A));
522 H->Setup(M, 0, H->GetNumLevels());
527 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
528 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
530 FilterMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Pattern,
Scalar dropTol)
const {
531 typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
537 RCP<GraphBase> filteredGraph;
543 level.
Set<RCP<Matrix> >(
"A", rcpFromRef(Pattern));
545 RCP<AmalgamationFactory> amalgFactory = rcp(
new AmalgamationFactory());
547 RCP<CoalesceDropFactory> dropFactory = rcp(
new CoalesceDropFactory());
548 ParameterList dropParams = *(dropFactory->GetValidParameterList());
549 dropParams.set(
"lightweight wrap",
true);
550 dropParams.set(
"aggregation: drop scheme",
"classical");
551 dropParams.set(
"aggregation: drop tol", dropTol);
553 dropFactory->SetParameterList(dropParams);
554 dropFactory->SetFactory(
"UnAmalgamationInfo", amalgFactory);
557 level.
Request(
"Graph", dropFactory.get());
558 dropFactory->Build(level);
560 level.
Get(
"Graph", filteredGraph, dropFactory.get());
563 RCP<Matrix> filteredA;
569 level.
Set(
"A", rcpFromRef(A));
570 level.
Set(
"Graph", filteredGraph);
571 level.
Set(
"Filtering",
true);
573 RCP<FilteredAFactory> filterFactory = rcp(
new FilteredAFactory());
574 ParameterList filterParams = *(filterFactory->GetValidParameterList());
577 filterParams.set(
"filtered matrix: reuse graph",
false);
578 filterParams.set(
"filtered matrix: use lumping",
false);
579 filterFactory->SetParameterList(filterParams);
582 level.
Request(
"A", filterFactory.get());
583 filterFactory->Build(level);
585 level.
Get(
"A", filteredA, filterFactory.get());
589 filteredA->resumeFill();
590 size_t numRows = filteredA->getRowMap()->getLocalNumElements();
591 for (
size_t i = 0; i < numRows; i++) {
592 ArrayView<const LO> inds;
593 ArrayView<const SC> vals;
594 filteredA->getLocalRowView(i, inds, vals);
596 size_t nnz = inds.size();
598 Array<SC> valsNew = vals;
601 SC diag = Teuchos::ScalarTraits<SC>::zero();
602 for (
size_t j = 0; j < nnz; j++) {
604 if (inds[j] == Teuchos::as<int>(i))
608 "No diagonal found");
612 valsNew[diagIndex] -= diag;
614 filteredA->replaceLocalValues(i, inds, valsNew);
616 filteredA->fillComplete();
621 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
625 typedef MueLu::BlockedPFactory <SC,LO,GO,NO> BlockedPFactory;
626 typedef MueLu::GenericRFactory <SC,LO,GO,NO> GenericRFactory;
627 typedef MueLu::BlockedRAPFactory <SC,LO,GO,NO> BlockedRAPFactory;
628 typedef MueLu::SmootherFactory <SC,LO,GO,NO> SmootherFactory;
629 typedef MueLu::BlockedDirectSolver <SC,LO,GO,NO> BlockedDirectSolver;
630 typedef MueLu::FactoryManager <SC,LO,GO,NO> FactoryManager;
635 RCP<FactoryManager> M11 = rcp(
new FactoryManager()), M22 = rcp(
new FactoryManager());
636 M11->SetKokkosRefactor(paramList.get<
bool>(
"use kokkos refactor"));
637 M22->SetKokkosRefactor(paramList.get<
bool>(
"use kokkos refactor"));
638 SetBlockDependencyTree(*M11, 0, 0,
"velocity", paramList);
639 SetBlockDependencyTree(*M22, 1, 1,
"pressure", paramList);
641 RCP<BlockedPFactory> PFact = rcp(
new BlockedPFactory());
642 ParameterList pParamList = *(PFact->GetValidParameterList());
643 pParamList.set(
"backwards",
true);
644 PFact->SetParameterList(pParamList);
645 PFact->AddFactoryManager(M11);
646 PFact->AddFactoryManager(M22);
649 RCP<GenericRFactory> RFact = rcp(
new GenericRFactory());
650 RFact->SetFactory(
"P", PFact);
653 RCP<MueLu::Factory > AcFact = rcp(
new BlockedRAPFactory());
654 AcFact->SetFactory(
"R", RFact);
655 AcFact->SetFactory(
"P", PFact);
661 RCP<MueLu::Factory> coarseFact = rcp(
new SmootherFactory(rcp(
new BlockedDirectSolver()), Teuchos::null));
666 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
670 typedef MueLu::ConstraintFactory <SC,LO,GO,NO> ConstraintFactory;
671 typedef MueLu::EminPFactory <SC,LO,GO,NO> EminPFactory;
672 typedef MueLu::GenericRFactory <SC,LO,GO,NO> GenericRFactory;
673 typedef MueLu::PatternFactory <SC,LO,GO,NO> PatternFactory;
674 typedef MueLu::Q2Q1PFactory <SC,LO,GO,NO> Q2Q1PFactory;
675 typedef MueLu::Q2Q1uPFactory <SC,LO,GO,NO> Q2Q1uPFactory;
676 typedef MueLu::SubBlockAFactory <SC,LO,GO,NO> SubBlockAFactory;
678 RCP<SubBlockAFactory> AFact = rcp(
new SubBlockAFactory());
680 AFact->SetParameter(
"block row", Teuchos::ParameterEntry(row));
681 AFact->SetParameter(
"block col", Teuchos::ParameterEntry(col));
684 RCP<MueLu::Factory> Q2Q1Fact;
686 const bool isStructured =
false;
689 Q2Q1Fact = rcp(
new Q2Q1PFactory);
692 Q2Q1Fact = rcp(
new Q2Q1uPFactory);
693 ParameterList q2q1ParamList = *(Q2Q1Fact->GetValidParameterList());
694 q2q1ParamList.set(
"mode", mode);
695 if (paramList.isParameter(
"dump status"))
696 q2q1ParamList.set(
"dump status", paramList.get<
bool>(
"dump status"));
697 if (paramList.isParameter(
"phase2"))
698 q2q1ParamList.set(
"phase2", paramList.get<
bool>(
"phase2"));
699 if (paramList.isParameter(
"tau_2"))
700 q2q1ParamList.set(
"tau_2", paramList.get<
double>(
"tau_2"));
701 Q2Q1Fact->SetParameterList(q2q1ParamList);
703 Q2Q1Fact->SetFactory(
"A", AFact);
706 RCP<PatternFactory> patternFact = rcp(
new PatternFactory);
707 ParameterList patternParams = *(patternFact->GetValidParameterList());
710 patternParams.set(
"emin: pattern order", 0);
711 patternFact->SetParameterList(patternParams);
712 patternFact->SetFactory(
"A", AFact);
713 patternFact->SetFactory(
"P", Q2Q1Fact);
716 RCP<ConstraintFactory> CFact = rcp(
new ConstraintFactory);
717 CFact->SetFactory(
"Ppattern", patternFact);
720 RCP<EminPFactory> EminPFact = rcp(
new EminPFactory());
721 ParameterList eminParams = *(EminPFact->GetValidParameterList());
722 if (paramList.isParameter(
"emin: num iterations"))
723 eminParams.set(
"emin: num iterations", paramList.get<
int>(
"emin: num iterations"));
724 if (mode ==
"pressure") {
725 eminParams.set(
"emin: iterative method",
"cg");
727 eminParams.set(
"emin: iterative method",
"gmres");
728 if (paramList.isParameter(
"emin: iterative method"))
729 eminParams.set(
"emin: iterative method", paramList.get<std::string>(
"emin: iterative method"));
731 EminPFact->SetParameterList(eminParams);
732 EminPFact->SetFactory(
"A", AFact);
733 EminPFact->SetFactory(
"Constraint", CFact);
734 EminPFact->SetFactory(
"P", Q2Q1Fact);
737 if (mode ==
"velocity" && (!paramList.isParameter(
"velocity: use transpose") || paramList.get<
bool>(
"velocity: use transpose") ==
false)) {
740 RCP<GenericRFactory> RFact = rcp(
new GenericRFactory());
741 RFact->SetFactory(
"P", EminPFact);
746 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
747 RCP<MueLu::FactoryBase>
749 GetSmoother(
const std::string& type,
const ParameterList& paramList,
bool coarseSolver)
const {
750 typedef Teuchos::ParameterEntry ParameterEntry;
752 typedef MueLu::BlockedDirectSolver <SC,LO,GO,NO> BlockedDirectSolver;
753 typedef MueLu::BraessSarazinSmoother <SC,LO,GO,NO> BraessSarazinSmoother;
754 typedef MueLu::DirectSolver <SC,LO,GO,NO> DirectSolver;
755 typedef MueLu::FactoryManager <SC,LO,GO,NO> FactoryManager;
757 typedef MueLu::SmootherFactory <SC,LO,GO,NO> SmootherFactory;
758 typedef MueLu::SmootherPrototype <SC,LO,GO,NO> SmootherPrototype;
759 typedef MueLu::TrilinosSmoother <SC,LO,GO,NO> TrilinosSmoother;
761 RCP<SmootherPrototype> smootherPrototype;
762 if (type ==
"none") {
763 return Teuchos::null;
765 }
else if (type ==
"vanka") {
767 ParameterList schwarzList;
768 schwarzList.set(
"schwarz: overlap level", as<int>(0));
769 schwarzList.set(
"schwarz: zero starting solution",
false);
770 schwarzList.set(
"subdomain solver name",
"Block_Relaxation");
772 ParameterList& innerSolverList = schwarzList.sublist(
"subdomain solver parameters");
773 innerSolverList.set(
"partitioner: type",
"user");
774 innerSolverList.set(
"partitioner: overlap",
MUELU_GPD(
"partitioner: overlap",
int, 1));
775 innerSolverList.set(
"relaxation: type",
MUELU_GPD(
"relaxation: type", std::string,
"Gauss-Seidel"));
776 innerSolverList.set(
"relaxation: sweeps",
MUELU_GPD(
"relaxation: sweeps",
int, 1));
777 innerSolverList.set(
"relaxation: damping factor",
MUELU_GPD(
"relaxation: damping factor",
double, 0.5));
778 innerSolverList.set(
"relaxation: zero starting solution",
false);
781 std::string ifpackType =
"SCHWARZ";
783 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, schwarzList));
785 }
else if (type ==
"schwarz") {
787 std::string ifpackType =
"SCHWARZ";
789 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, paramList));
791 }
else if (type ==
"braess-sarazin") {
794 bool lumping =
MUELU_GPD(
"bs: lumping",
bool,
false);
796 RCP<SchurComplementFactory> schurFact = rcp(
new SchurComplementFactory());
797 schurFact->SetParameter(
"omega", ParameterEntry(omega));
798 schurFact->SetParameter(
"lumping", ParameterEntry(lumping));
802 RCP<SmootherPrototype> schurSmootherPrototype;
803 std::string schurSmootherType = (paramList.isParameter(
"schur smoother: type") ? paramList.get<std::string>(
"schur smoother: type") :
"RELAXATION");
804 if (schurSmootherType ==
"RELAXATION") {
805 ParameterList schurSmootherParams = paramList.sublist(
"schur smoother: params");
807 schurSmootherPrototype = rcp(
new TrilinosSmoother(schurSmootherType, schurSmootherParams));
809 schurSmootherPrototype = rcp(
new DirectSolver());
811 schurSmootherPrototype->SetFactory(
"A", schurFact);
813 RCP<SmootherFactory> schurSmootherFact = rcp(
new SmootherFactory(schurSmootherPrototype));
816 RCP<FactoryManager> braessManager = rcp(
new FactoryManager());
817 braessManager->SetFactory(
"A", schurFact);
818 braessManager->SetFactory(
"Smoother", schurSmootherFact);
819 braessManager->SetFactory(
"PreSmoother", schurSmootherFact);
820 braessManager->SetFactory(
"PostSmoother", schurSmootherFact);
821 braessManager->SetIgnoreUserData(
true);
823 smootherPrototype = rcp(
new BraessSarazinSmoother());
824 smootherPrototype->SetParameter(
"Sweeps", ParameterEntry(
MUELU_GPD(
"bs: sweeps",
int, 1)));
825 smootherPrototype->SetParameter(
"lumping", ParameterEntry(lumping));
826 smootherPrototype->SetParameter(
"Damping factor", ParameterEntry(omega));
827 smootherPrototype->SetParameter(
"q2q1 mode", ParameterEntry(
true));
828 rcp_dynamic_cast<BraessSarazinSmoother>(smootherPrototype)->AddFactoryManager(braessManager, 0);
830 }
else if (type ==
"ilu") {
831 std::string ifpackType =
"RILUK";
833 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, paramList));
835 }
else if (type ==
"direct") {
836 smootherPrototype = rcp(
new BlockedDirectSolver());
842 return coarseSolver ? rcp(
new SmootherFactory(smootherPrototype, Teuchos::null)) : rcp(
new SmootherFactory(smootherPrototype));
845 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
846 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
848 typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
849 typedef Xpetra::CrsMatrixWrap<SC,LO,GO,NO> CrsMatrixWrap;
850 typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
852 const CrsMatrixWrap& Awrap =
dynamic_cast<const CrsMatrixWrap&
>(A);
854 ArrayRCP<const size_t> iaA;
855 ArrayRCP<const LO> jaA;
856 ArrayRCP<const SC> valA;
857 Awrap.getCrsMatrix()->getAllValues(iaA, jaA, valA);
859 ArrayRCP<size_t> iaB (iaA .size());
860 ArrayRCP<LO> jaB (jaA .size());
861 ArrayRCP<SC> valB(valA.size());
862 for (
int i = 0; i < iaA .size(); i++) iaB [i] = iaA[i];
863 for (
int i = 0; i < jaA .size(); i++) jaB [i] = jaA[i];
864 for (
int i = 0; i < valA.size(); i++) valB[i] = Teuchos::ScalarTraits<SC>::magnitude(valA[i]);
866 RCP<Matrix> B = rcp(
new CrsMatrixWrap(A.getRowMap(), A.getColMap(), 0));
867 RCP<CrsMatrix> Bcrs = rcp_dynamic_cast<CrsMatrixWrap>(B)->getCrsMatrix();
868 Bcrs->setAllValues(iaB, jaB, valB);
869 Bcrs->expertStaticFillComplete(A.getDomainMap(), A.getRangeMap());
875 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
877 return "Thyra::MueLuTpetraQ2Q1PreconditionerFactory";
Various adapters that will create a MueLu preconditioner that is a Tpetra::Operator.
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
#define MUELU_GPD(name, type, defaultValue)
AmalgamationFactory for subblocks of strided map based amalgamation data.
Factory for creating a graph based on a given matrix.
Exception throws to report errors in the internal logical of the program.
This class specifies the default factory that should generate some data on a Level if the data does n...
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.
Factory for building filtered matrices using filtered graphs.
MueLu representation of a graph.
Class that holds all level-specific information.
void SetLevelID(int levelID)
Set level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
static const RCP< const NoFactory > getRCP()
Static Get() functions.
Factory for building the Schur Complement for a 2x2 block matrix.
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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
MueLuTpetraQ2Q1PreconditionerFactory()
Teuchos::RCP< PreconditionerBase< SC > > createPrec() const
void SetDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, const ParameterList ¶mList) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::RCP< Xpetra::Matrix< SC, LO, GO, NO > > Absolute(const Xpetra::Matrix< SC, LO, GO, NO > &A) const
void uninitializePrec(PreconditionerBase< SC > *prec, Teuchos::RCP< const LinearOpSourceBase< SC > > *fwdOp, ESupportSolveUse *supportSolveUse) const
std::string description() const
void SetBlockDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, LO row, LO col, const std::string &mode, const ParameterList ¶mList) const
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
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
RCP< MueLu::FactoryBase > GetSmoother(const std::string &type, const ParameterList ¶mList, bool coarseSolver) const
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< SC > > &fwdOp, PreconditionerBase< SC > *prec, const ESupportSolveUse supportSolveUse) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
bool isCompatible(const LinearOpSourceBase< SC > &fwdOp) const
RCP< Xpetra::Matrix< SC, LO, GO, NO > > TpetraCrs_To_XpetraMatrix(const Teuchos::RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > &Atpetra)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.