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;
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()->getNodeNumElements();
315 Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getNodeNumElements();
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->getNodeNumElements();
322 Xpetra::global_size_t numCols2 = domainMap2->getNodeNumElements();
323 ArrayView<const GO> rangeElem2 = rangeMap2->getNodeElementList();
324 ArrayView<const GO> domainElem2 = domainMap2->getNodeElementList();
325 ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getNodeElementList();
326 ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getNodeElementList();
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->getNodeMaxNumRowEntries());
345 RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getNodeMaxNumRowEntries());
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()->getNodeNumElements()); ++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()->getNodeElementList();
437 Array<GO> presElem = A_21_crs->getRangeMap()->getNodeElementList();
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()->getNodeNumElements();
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>
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>
674 typedef MueLu::Q2Q1PFactory <SC,LO,GO,NO> Q2Q1PFactory;
675 typedef MueLu::Q2Q1uPFactory <SC,LO,GO,NO> Q2Q1uPFactory;
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;
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";
883 #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
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
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()
MueLu::DefaultScalar Scalar
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.