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;
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 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::MapExtractor <SC,LO,GO,NO> MapExtractor;
285 typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
286 typedef Xpetra::MapFactory <LO,GO,NO> MapFactory;
287 typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
288 typedef Xpetra::MatrixFactory <SC,LO,GO,NO> MatrixFactory;
289 typedef Xpetra::StridedMapFactory <LO,GO,NO> StridedMapFactory;
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);
316 Xpetra::global_size_t numVel = A_11->getRowMap()->getNodeNumElements();
317 Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getNodeNumElements();
321 RCP<const Map> domainMap2 = tmp_A_12->getDomainMap();
322 RCP<const Map> rangeMap2 = tmp_A_21->getRangeMap();
323 Xpetra::global_size_t numRows2 = rangeMap2->getNodeNumElements();
324 Xpetra::global_size_t numCols2 = domainMap2->getNodeNumElements();
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();
330 Xpetra::UnderlyingLib lib = domainMap2->lib();
331 GO indexBase = domainMap2->getIndexBase();
333 Array<GO> newRowElem2(numRows2, 0);
334 for (Xpetra::global_size_t i = 0; i < numRows2; i++)
335 newRowElem2[i] = numVel + rangeElem2[i];
337 RCP<const Map> newRangeMap2 = MapFactory::Build(lib, numRows2, newRowElem2, indexBase, comm);
340 Array<GO> newColElem2(numCols2, 0);
341 for (Xpetra::global_size_t i = 0; i < numCols2; i++)
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);
421 RCP<Matrix> BBt = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21,
false, *A_12,
false, out);
422 RCP<Matrix> BBt_abs = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21_abs,
false, *A_12_abs,
false, out);
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();
441 for (Xpetra::global_size_t i = 0 ; i < numVel; i++) elementList[i] = velElem[i];
442 for (Xpetra::global_size_t i = numVel; i < numVel+numPres; i++) elementList[i] = presElem[i-numVel];
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);
468 finestLevel->Set(
"CoordinatesVelocity", Xpetra::toXpetra(velCoords));
469 finestLevel->Set(
"CoordinatesPressure", Xpetra::toXpetra(presCoords));
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);
491 Xpetra::IO<SC,LO,GO,NO>::Write(
"Pp_l" +
MueLu::toString(i) +
".mm", *Pp);
492 Xpetra::IO<SC,LO,GO,NO>::Write(
"Pv_l" +
MueLu::toString(i) +
".mm", *Pv);
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> >
532 FilterMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Pattern, Scalar dropTol)
const {
533 typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
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> >
847 typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
848 typedef Xpetra::CrsMatrixWrap<SC,LO,GO,NO> CrsMatrixWrap;
849 typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
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]);
865 RCP<Matrix> B = rcp(
new CrsMatrixWrap(A.getRowMap(), A.getColMap(), 0, Xpetra::StaticProfile));
866 RCP<CrsMatrix> Bcrs = rcp_dynamic_cast<CrsMatrixWrap>(B)->getCrsMatrix();
867 Bcrs->setAllValues(iaB, jaB, valB);
868 Bcrs->expertStaticFillComplete(A.getDomainMap(), A.getRangeMap());
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.
bool isCompatible(const LinearOpSourceBase< SC > &fwdOp) const
Class that encapsulates external library smoothers.
Factory for building coarse matrices.
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.
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< 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.
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
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.