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" 107 #define MUELU_GPD(name, type, defaultValue) \ 108 (paramList.isParameter(name) ? paramList.get<type>(name) : defaultValue) 112 using Teuchos::rcp_const_cast;
113 using Teuchos::rcp_dynamic_cast;
114 using Teuchos::ParameterList;
115 using Teuchos::ArrayView;
116 using Teuchos::ArrayRCP;
118 using Teuchos::Array;
121 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
128 typedef Thyra ::TpetraLinearOp<SC,LO,GO,NO> ThyraTpetraLinOp;
129 typedef Tpetra::Operator <SC,LO,GO,NO> TpetraLinOp;
130 typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TpetraCrsMat;
132 const RCP<const LinearOpBase<SC> > fwdOp = fwdOpSrc.getOp();
133 const RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = rcp_dynamic_cast<
const ThyraTpetraLinOp>(fwdOp);
134 const RCP<const TpetraLinOp> tpetraFwdOp = Teuchos::nonnull(thyraTpetraFwdOp) ? thyraTpetraFwdOp->getConstTpetraOperator() : Teuchos::null;
135 const RCP<const TpetraCrsMat> tpetraFwdCrsMat = rcp_dynamic_cast<
const TpetraCrsMat>(tpetraFwdOp);
137 return Teuchos::nonnull(tpetraFwdCrsMat);
140 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
141 RCP<PreconditionerBase<Scalar> >
143 return rcp(
new DefaultPreconditioner<SC>);
146 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
148 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> > &fwdOpSrc, PreconditionerBase<Scalar> *prec,
const ESupportSolveUse supportSolveUse)
const {
150 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
151 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
152 TEUCHOS_ASSERT(prec);
155 const RCP<const LinearOpBase<SC> > fwdOp = fwdOpSrc->getOp();
156 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
158 typedef Thyra::TpetraLinearOp<SC,LO,GO,NO> ThyraTpetraLinOp;
159 const RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = rcp_dynamic_cast<
const ThyraTpetraLinOp>(fwdOp);
160 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraTpetraFwdOp));
162 typedef Tpetra::Operator<SC,LO,GO,NO> TpetraLinOp;
163 const RCP<const TpetraLinOp> tpetraFwdOp = thyraTpetraFwdOp->getConstTpetraOperator();
164 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdOp));
166 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraCrsMat;
167 const RCP<const TpetraCrsMat> tpetraFwdCrsMat = rcp_dynamic_cast<
const TpetraCrsMat>(tpetraFwdOp);
168 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdCrsMat));
171 const Teuchos::Ptr<DefaultPreconditioner<SC> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<SC> *
>(prec));
172 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
175 const RCP<TpetraCrsMat> tpetraFwdCrsMatNonConst = rcp_const_cast<TpetraCrsMat>(tpetraFwdCrsMat);
180 ParameterList paramList = *paramList_;
182 typedef Tpetra::MultiVector<SC,LO,GO,NO> MultiVector;
183 RCP<MultiVector> coords, nullspace, velCoords, presCoords;
185 Teko::LinearOp thA11, thA12, thA21, thA11_9Pt;
186 if (paramList.isType<RCP<MultiVector> >(
"Coordinates")) { coords = paramList.get<RCP<MultiVector> >(
"Coordinates"); paramList.remove(
"Coordinates"); }
187 if (paramList.isType<RCP<MultiVector> >(
"Nullspace")) { nullspace = paramList.get<RCP<MultiVector> >(
"Nullspace"); paramList.remove(
"Nullspace"); }
188 if (paramList.isType<RCP<MultiVector> >(
"Velcoords")) { velCoords = paramList.get<RCP<MultiVector> >(
"Velcoords"); paramList.remove(
"Velcoords"); }
189 if (paramList.isType<RCP<MultiVector> >(
"Prescoords")) { presCoords = paramList.get<RCP<MultiVector> >(
"Prescoords"); paramList.remove(
"Prescoords"); }
190 if (paramList.isType<ArrayRCP<LO> > (
"p2vMap")) { p2vMap = paramList.get<ArrayRCP<LO> > (
"p2vMap"); paramList.remove(
"p2vMap"); }
191 if (paramList.isType<Teko::LinearOp> (
"A11")) { thA11 = paramList.get<Teko::LinearOp> (
"A11"); paramList.remove(
"A11"); }
192 if (paramList.isType<Teko::LinearOp> (
"A12")) { thA12 = paramList.get<Teko::LinearOp> (
"A12"); paramList.remove(
"A12"); }
193 if (paramList.isType<Teko::LinearOp> (
"A21")) { thA21 = paramList.get<Teko::LinearOp> (
"A21"); paramList.remove(
"A21"); }
194 if (paramList.isType<Teko::LinearOp> (
"A11_9Pt")) { thA11_9Pt = paramList.get<Teko::LinearOp> (
"A11_9Pt"); paramList.remove(
"A11_9Pt"); }
197 const RCP<MueLuOperator> mueluPrecOp = Q2Q1MkPrecond(paramList, velCoords, presCoords, p2vMap, thA11, thA12, thA21, thA11_9Pt);
199 const RCP<LinearOpBase<SC> > thyraPrecOp = Thyra::createLinearOp(RCP<TpetraLinOp>(mueluPrecOp));
200 defaultPrec->initializeUnspecified(thyraPrecOp);
203 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
205 uninitializePrec(PreconditionerBase<Scalar> *prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
207 TEUCHOS_ASSERT(prec);
210 const Teuchos::Ptr<DefaultPreconditioner<SC> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<SC> *
>(prec));
211 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
215 *fwdOp = Teuchos::null;
218 if (supportSolveUse) {
220 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
223 defaultPrec->uninitialize();
228 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
230 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
231 paramList_ = paramList;
234 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
241 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
244 RCP<ParameterList> savedParamList = paramList_;
245 paramList_ = Teuchos::null;
246 return savedParamList;
249 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
250 RCP<const ParameterList>
255 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
256 RCP<const ParameterList>
258 static RCP<const ParameterList> validPL;
260 if (validPL.is_null())
261 validPL = rcp(
new ParameterList());
266 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
267 RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
270 const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& velCoords,
271 const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& presCoords,
272 const ArrayRCP<LocalOrdinal>& p2vMap,
273 const Teko::LinearOp& thA11,
const Teko::LinearOp& thA12,
const Teko::LinearOp& thA21,
const Teko::LinearOp& thA11_9Pt)
const 277 typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TP_Crs;
278 typedef Tpetra::Operator <SC,LO,GO,NO> TP_Op;
280 typedef Xpetra::BlockedCrsMatrix <SC,LO,GO,NO> BlockedCrsMatrix;
281 typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
282 typedef Xpetra::CrsMatrixWrap <SC,LO,GO,NO> CrsMatrixWrap;
283 typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
284 typedef Xpetra::MapExtractor <SC,LO,GO,NO> MapExtractor;
285 typedef Xpetra::Map <LO,GO,NO> Map;
286 typedef Xpetra::MapExtractor <SC,LO,GO,NO> MapExtractor;
287 typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
288 typedef Xpetra::MapFactory <LO,GO,NO> MapFactory;
289 typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
290 typedef Xpetra::MatrixFactory <SC,LO,GO,NO> MatrixFactory;
291 typedef Xpetra::StridedMapFactory <LO,GO,NO> StridedMapFactory;
295 const RCP<const Teuchos::Comm<int> > comm = velCoords->getMap()->getComm();
298 RCP<Thyra::LinearOpBase<SC> > ThNonConstA11 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11);
299 RCP<Thyra::LinearOpBase<SC> > ThNonConstA21 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA21);
300 RCP<Thyra::LinearOpBase<SC> > ThNonConstA12 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA12);
301 RCP<Thyra::LinearOpBase<SC> > ThNonConstA11_9Pt = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11_9Pt);
303 RCP<TP_Op> TpetA11 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11);
304 RCP<TP_Op> TpetA21 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA21);
305 RCP<TP_Op> TpetA12 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA12);
306 RCP<TP_Op> TpetA11_9Pt = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11_9Pt);
308 RCP<TP_Crs> TpetCrsA11 = rcp_dynamic_cast<TP_Crs>(TpetA11);
309 RCP<TP_Crs> TpetCrsA21 = rcp_dynamic_cast<TP_Crs>(TpetA21);
310 RCP<TP_Crs> TpetCrsA12 = rcp_dynamic_cast<TP_Crs>(TpetA12);
311 RCP<TP_Crs> TpetCrsA11_9Pt = rcp_dynamic_cast<TP_Crs>(TpetA11_9Pt);
318 Xpetra::global_size_t numVel = A_11->getRowMap()->getNodeNumElements();
319 Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getNodeNumElements();
323 RCP<const Map> domainMap2 = tmp_A_12->getDomainMap();
324 RCP<const Map> rangeMap2 = tmp_A_21->getRangeMap();
325 Xpetra::global_size_t numRows2 = rangeMap2->getNodeNumElements();
326 Xpetra::global_size_t numCols2 = domainMap2->getNodeNumElements();
327 ArrayView<const GO> rangeElem2 = rangeMap2->getNodeElementList();
328 ArrayView<const GO> domainElem2 = domainMap2->getNodeElementList();
329 ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getNodeElementList();
330 ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getNodeElementList();
332 Xpetra::UnderlyingLib lib = domainMap2->lib();
333 GO indexBase = domainMap2->getIndexBase();
335 Array<GO> newRowElem2(numRows2, 0);
336 for (Xpetra::global_size_t i = 0; i < numRows2; i++)
337 newRowElem2[i] = numVel + rangeElem2[i];
339 RCP<const Map> newRangeMap2 = MapFactory::Build(lib, numRows2, newRowElem2, indexBase, comm);
342 Array<GO> newColElem2(numCols2, 0);
343 for (Xpetra::global_size_t i = 0; i < numCols2; i++)
344 newColElem2[i] = numVel + domainElem2[i];
346 RCP<const Map> newDomainMap2 = MapFactory::Build(lib, numCols2, newColElem2, indexBase, comm);
348 RCP<Matrix> A_12 = MatrixFactory::Build(tmp_A_12->getRangeMap(), newDomainMap2, tmp_A_12->getNodeMaxNumRowEntries());
349 RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getNodeMaxNumRowEntries());
351 RCP<CrsMatrix> A_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_11) ->getCrsMatrix();
352 RCP<CrsMatrix> A_12_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_12) ->getCrsMatrix();
353 RCP<CrsMatrix> A_21_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_21) ->getCrsMatrix();
354 RCP<CrsMatrix> A_11_crs_9Pt = rcp_dynamic_cast<CrsMatrixWrap>(A_11_9Pt)->getCrsMatrix();
357 RCP<Matrix> A_22 = MatrixFactory::Build(newRangeMap2, newDomainMap2, 1);
358 RCP<CrsMatrix> A_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_22) ->getCrsMatrix();
361 Array<SC> smallVal(1, 1.0e-9);
366 ArrayView<const LO> inds;
367 ArrayView<const SC> vals;
368 for (
LO row = 0; row < as<LO>(numRows2); ++row) {
369 tmp_A_21->getLocalRowView(row, inds, vals);
371 size_t nnz = inds.size();
372 Array<GO> newInds(nnz, 0);
373 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
374 newInds[colID] = colElem1[inds[colID]];
376 A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
377 A_22_crs->insertGlobalValues(newRowElem2[row], Array<LO>(1, newRowElem2[row]), smallVal);
379 A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
380 A_22_crs->fillComplete(newDomainMap2, newRangeMap2);
382 RCP<Matrix> A_22 = Teuchos::null;
383 RCP<CrsMatrix> A_22_crs = Teuchos::null;
385 ArrayView<const LO> inds;
386 ArrayView<const SC> vals;
387 for (
LO row = 0; row < as<LO>(numRows2); ++row) {
388 tmp_A_21->getLocalRowView(row, inds, vals);
390 size_t nnz = inds.size();
391 Array<GO> newInds(nnz, 0);
392 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
393 newInds[colID] = colElem1[inds[colID]];
395 A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
397 A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
402 for (
LO row = 0; row < as<LO>(tmp_A_12->getRowMap()->getNodeNumElements()); ++row) {
403 tmp_A_12->getLocalRowView(row, inds, vals);
405 size_t nnz = inds.size();
406 Array<GO> newInds(nnz, 0);
407 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
408 newInds[colID] = newColElem2[inds[colID]];
410 A_12_crs->insertGlobalValues(rowElem1[row], newInds, vals);
412 A_12_crs->fillComplete(newDomainMap2, tmp_A_12->getRangeMap());
414 RCP<Matrix> A_12_abs = Absolute(*A_12);
415 RCP<Matrix> A_21_abs = Absolute(*A_21);
420 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
421 Teuchos::FancyOStream& out = *fancy;
422 out.setOutputToRootOnly(0);
423 RCP<Matrix> BBt = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21,
false, *A_12,
false, out);
424 RCP<Matrix> BBt_abs = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21_abs,
false, *A_12_abs,
false, out);
426 SC dropTol = (paramList.get<
int>(
"useFilters") ? 0.06 : 0.00);
427 RCP<Matrix> filteredA = FilterMatrix(*A_11, *A_11, dropTol);
428 RCP<Matrix> filteredB = FilterMatrix(*BBt, *BBt_abs, dropTol);
430 RCP<CrsMatrix> fA_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredA)->getCrsMatrix();
431 RCP<CrsMatrix> fA_12_crs = Teuchos::null;
432 RCP<CrsMatrix> fA_21_crs = Teuchos::null;
433 RCP<CrsMatrix> fA_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredB)->getCrsMatrix();
436 std::vector<size_t> stridingInfo(1, 1);
437 int stridedBlockId = -1;
439 Array<GO> elementList(numVel+numPres);
440 Array<GO> velElem = A_12_crs->getRangeMap()->getNodeElementList();
441 Array<GO> presElem = A_21_crs->getRangeMap()->getNodeElementList();
443 for (Xpetra::global_size_t i = 0 ; i < numVel; i++) elementList[i] = velElem[i];
444 for (Xpetra::global_size_t i = numVel; i < numVel+numPres; i++) elementList[i] = presElem[i-numVel];
445 RCP<const Map> fullMap = StridedMapFactory::Build(Xpetra::UseTpetra, numVel+numPres, elementList(), indexBase, stridingInfo, comm);
447 std::vector<RCP<const Map> > partMaps(2);
448 partMaps[0] = StridedMapFactory::Build(Xpetra::UseTpetra, numVel, velElem, indexBase, stridingInfo, comm);
449 partMaps[1] = StridedMapFactory::Build(Xpetra::UseTpetra, numPres, presElem, indexBase, stridingInfo, comm, stridedBlockId, numVel);
452 RCP<const MapExtractor> mapExtractor = MapExtractorFactory::Build(fullMap, partMaps);
453 RCP<BlockedCrsMatrix> fA = rcp(
new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
454 fA->setMatrix(0, 0, fA_11_crs);
455 fA->setMatrix(0, 1, fA_12_crs);
456 fA->setMatrix(1, 0, fA_21_crs);
457 fA->setMatrix(1, 1, fA_22_crs);
464 SetDependencyTree(M, paramList);
466 RCP<Hierarchy> H = rcp(
new Hierarchy);
467 RCP<MueLu::Level> finestLevel = H->GetLevel(0);
468 finestLevel->Set(
"A", rcp_dynamic_cast<Matrix>(fA));
469 finestLevel->Set(
"p2vMap", p2vMap);
470 finestLevel->Set(
"CoordinatesVelocity", Xpetra::toXpetra(velCoords));
471 finestLevel->Set(
"CoordinatesPressure", Xpetra::toXpetra(presCoords));
472 finestLevel->Set(
"AForPat", A_11_9Pt);
473 H->SetMaxCoarseSize(
MUELU_GPD(
"coarse: max size",
int, 1));
475 #ifdef IMPLICIT_TRANSPOSE 476 out <<
"Using implicit transpose" << std::endl;
478 H->SetImplicitTranspose(
true);
489 H->Keep(
"Ptent", M.
GetFactory(
"Ptent").get());
490 H->Setup(M, 0,
MUELU_GPD(
"max levels",
int, 3));
493 for (
int i = 1; i < H->GetNumLevels(); i++) {
494 RCP<Matrix> P = H->GetLevel(i)->template Get<RCP<Matrix> >(
"P");
495 RCP<BlockedCrsMatrix> Pcrs = rcp_dynamic_cast<BlockedCrsMatrix>(P);
496 RCP<CrsMatrix> Ppcrs = Pcrs->getMatrix(1,1);
497 RCP<Matrix> Pp = rcp(
new CrsMatrixWrap(Ppcrs));
498 RCP<CrsMatrix> Pvcrs = Pcrs->getMatrix(0,0);
499 RCP<Matrix> Pv = rcp(
new CrsMatrixWrap(Pvcrs));
501 Xpetra::IO<SC,LO,GO,NO>::Write(
"Pp_l" +
MueLu::toString(i) +
".mm", *Pp);
502 Xpetra::IO<SC,LO,GO,NO>::Write(
"Pv_l" +
MueLu::toString(i) +
".mm", *Pv);
509 std::string smootherType =
MUELU_GPD(
"smoother: type", std::string,
"vanka");
510 ParameterList smootherParams;
511 if (paramList.isSublist(
"smoother: params"))
512 smootherParams = paramList.sublist(
"smoother: params");
513 M.
SetFactory(
"Smoother", GetSmoother(smootherType, smootherParams,
false));
515 std::string coarseType =
MUELU_GPD(
"coarse: type", std::string,
"direct");
516 ParameterList coarseParams;
517 if (paramList.isSublist(
"coarse: params"))
518 coarseParams = paramList.sublist(
"coarse: params");
519 M.
SetFactory(
"CoarseSolver", GetSmoother(coarseType, coarseParams,
true));
521 #ifdef HAVE_MUELU_DEBUG 525 RCP<BlockedCrsMatrix> A = rcp(
new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
526 A->setMatrix(0, 0, A_11_crs);
527 A->setMatrix(0, 1, A_12_crs);
528 A->setMatrix(1, 0, A_21_crs);
529 A->setMatrix(1, 1, A_22_crs);
532 H->GetLevel(0)->Set(
"A", rcp_dynamic_cast<Matrix>(A));
534 H->Setup(M, 0, H->GetNumLevels());
539 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
540 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
542 FilterMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Pattern, Scalar dropTol)
const {
543 typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
550 RCP<GraphBase> filteredGraph;
556 level.
Set<RCP<Matrix> >(
"A", rcpFromRef(Pattern));
558 RCP<AmalgamationFactory> amalgFactory = rcp(
new AmalgamationFactory());
560 RCP<CoalesceDropFactory> dropFactory = rcp(
new CoalesceDropFactory());
561 ParameterList dropParams = *(dropFactory->GetValidParameterList());
562 dropParams.set(
"lightweight wrap",
true);
563 dropParams.set(
"aggregation: drop scheme",
"classical");
564 dropParams.set(
"aggregation: drop tol", dropTol);
566 dropFactory->SetParameterList(dropParams);
567 dropFactory->SetFactory(
"UnAmalgamationInfo", amalgFactory);
570 level.
Request(
"Graph", dropFactory.get());
571 dropFactory->Build(level);
573 level.
Get(
"Graph", filteredGraph, dropFactory.get());
576 RCP<Matrix> filteredA;
582 level.
Set(
"A", rcpFromRef(A));
583 level.
Set(
"Graph", filteredGraph);
584 level.
Set(
"Filtering",
true);
586 RCP<FilteredAFactory> filterFactory = rcp(
new FilteredAFactory());
587 ParameterList filterParams = *(filterFactory->GetValidParameterList());
590 filterParams.set(
"filtered matrix: reuse graph",
false);
591 filterParams.set(
"filtered matrix: use lumping",
false);
592 filterFactory->SetParameterList(filterParams);
595 level.
Request(
"A", filterFactory.get());
596 filterFactory->Build(level);
598 level.
Get(
"A", filteredA, filterFactory.get());
602 filteredA->resumeFill();
603 size_t numRows = filteredA->getRowMap()->getNodeNumElements();
604 for (
size_t i = 0; i < numRows; i++) {
605 ArrayView<const LO> inds;
606 ArrayView<const SC> vals;
607 filteredA->getLocalRowView(i, inds, vals);
609 size_t nnz = inds.size();
611 Array<SC> valsNew = vals;
614 SC diag = Teuchos::ScalarTraits<SC>::zero();
615 for (
size_t j = 0; j < inds.size(); j++) {
621 "No diagonal found");
625 valsNew[diagIndex] -= diag;
627 filteredA->replaceLocalValues(i, inds, valsNew);
629 filteredA->fillComplete();
634 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
648 RCP<FactoryManager> M11 = rcp(
new FactoryManager()), M22 = rcp(
new FactoryManager());
649 SetBlockDependencyTree(*M11, 0, 0,
"velocity", paramList);
650 SetBlockDependencyTree(*M22, 1, 1,
"pressure", paramList);
652 RCP<BlockedPFactory> PFact = rcp(
new BlockedPFactory());
653 ParameterList pParamList = *(PFact->GetValidParameterList());
654 pParamList.set(
"backwards",
true);
655 PFact->SetParameterList(pParamList);
656 PFact->AddFactoryManager(M11);
657 PFact->AddFactoryManager(M22);
660 RCP<MueLu::Factory > AcFact = rcp(
new BlockedRAPFactory());
661 #ifdef IMPLICIT_TRANSPOSE 664 ParameterList RAPparams;
665 RAPparams.set(
"transpose: use implicit",
true);
666 AcFact->SetParameterList(RAPparams);
668 RCP<GenericRFactory> RFact = rcp(
new GenericRFactory());
669 RFact->SetFactory(
"P", PFact);
672 AcFact->SetFactory(
"R", RFact);
674 AcFact->SetFactory(
"P", PFact);
680 RCP<MueLu::Factory> coarseFact = rcp(
new SmootherFactory(rcp(
new BlockedDirectSolver()), Teuchos::null));
685 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
692 typedef MueLu::Q2Q1PFactory <SC,LO,GO,NO> Q2Q1PFactory;
693 typedef MueLu::Q2Q1uPFactory <SC,LO,GO,NO> Q2Q1uPFactory;
696 RCP<SubBlockAFactory> AFact = rcp(
new SubBlockAFactory());
698 AFact->SetParameter(
"block row", Teuchos::ParameterEntry(row));
699 AFact->SetParameter(
"block col", Teuchos::ParameterEntry(col));
702 RCP<MueLu::Factory> Q2Q1Fact;
704 const bool isStructured =
false;
707 Q2Q1Fact = rcp(
new Q2Q1PFactory);
710 Q2Q1Fact = rcp(
new Q2Q1uPFactory);
711 ParameterList q2q1ParamList = *(Q2Q1Fact->GetValidParameterList());
712 q2q1ParamList.set(
"mode", mode);
713 if (paramList.isParameter(
"dump status"))
714 q2q1ParamList.set(
"dump status", paramList.get<
bool>(
"dump status"));
715 if (paramList.isParameter(
"phase2"))
716 q2q1ParamList.set(
"phase2", paramList.get<
bool>(
"phase2"));
717 Q2Q1Fact->SetParameterList(q2q1ParamList);
719 Q2Q1Fact->SetFactory(
"A", AFact);
722 RCP<PatternFactory> patternFact = rcp(
new PatternFactory);
723 ParameterList patternParams = *(patternFact->GetValidParameterList());
726 patternParams.set(
"emin: pattern order", 0);
727 patternFact->SetParameterList(patternParams);
728 patternFact->SetFactory(
"A", AFact);
729 patternFact->SetFactory(
"P", Q2Q1Fact);
732 RCP<ConstraintFactory> CFact = rcp(
new ConstraintFactory);
733 CFact->SetFactory(
"Ppattern", patternFact);
736 RCP<EminPFactory> EminPFact = rcp(
new EminPFactory());
737 ParameterList eminParams = *(EminPFact->GetValidParameterList());
738 if (paramList.isParameter(
"emin: num iterations"))
739 eminParams.set(
"emin: num iterations", paramList.get<
int>(
"emin: num iterations"));
740 EminPFact->SetParameterList(eminParams);
741 EminPFact->SetFactory(
"A", AFact);
742 EminPFact->SetFactory(
"Constraint", CFact);
743 EminPFact->SetFactory(
"P", Q2Q1Fact);
747 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
748 RCP<MueLu::FactoryBase>
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 ==
"direct") {
786 smootherPrototype = rcp(
new BlockedDirectSolver());
788 }
else if (type ==
"braess-sarazin") {
791 bool lumping =
MUELU_GPD(
"bs: lumping",
bool,
false);
793 RCP<SchurComplementFactory> schurFact = rcp(
new SchurComplementFactory());
794 schurFact->SetParameter(
"omega", ParameterEntry(omega));
795 schurFact->SetParameter(
"lumping", ParameterEntry(lumping));
799 RCP<SmootherPrototype> schurSmootherPrototype;
800 std::string schurSmootherType = (paramList.isParameter(
"schur smoother: type") ? paramList.get<std::string>(
"schur smoother: type") :
"RELAXATION");
801 if (schurSmootherType ==
"RELAXATION") {
802 ParameterList schurSmootherParams = paramList.sublist(
"schur smoother: params");
804 schurSmootherPrototype = rcp(
new TrilinosSmoother(schurSmootherType, schurSmootherParams));
806 schurSmootherPrototype = rcp(
new DirectSolver());
808 schurSmootherPrototype->SetFactory(
"A", schurFact);
810 RCP<SmootherFactory> schurSmootherFact = rcp(
new SmootherFactory(schurSmootherPrototype));
813 RCP<FactoryManager> braessManager = rcp(
new FactoryManager());
814 braessManager->SetFactory(
"A", schurFact);
815 braessManager->SetFactory(
"Smoother", schurSmootherFact);
816 braessManager->SetFactory(
"PreSmoother", schurSmootherFact);
817 braessManager->SetFactory(
"PostSmoother", schurSmootherFact);
818 braessManager->SetIgnoreUserData(
true);
820 smootherPrototype = rcp(
new BraessSarazinSmoother());
821 smootherPrototype->SetParameter(
"Sweeps", ParameterEntry(
MUELU_GPD(
"bs: sweeps",
int, 1)));
822 smootherPrototype->SetParameter(
"lumping", ParameterEntry(lumping));
823 smootherPrototype->SetParameter(
"Damping factor", ParameterEntry(omega));
824 rcp_dynamic_cast<BraessSarazinSmoother>(smootherPrototype)->AddFactoryManager(braessManager, 0);
827 return coarseSolver ? rcp(
new SmootherFactory(smootherPrototype, Teuchos::null)) : rcp(
new SmootherFactory(smootherPrototype));
830 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
831 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
833 typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
834 typedef Xpetra::CrsMatrixWrap<SC,LO,GO,NO> CrsMatrixWrap;
835 typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
837 const CrsMatrixWrap& Awrap =
dynamic_cast<const CrsMatrixWrap&
>(A);
839 ArrayRCP<const size_t> iaA;
840 ArrayRCP<const LO> jaA;
841 ArrayRCP<const SC> valA;
842 Awrap.getCrsMatrix()->getAllValues(iaA, jaA, valA);
844 ArrayRCP<size_t> iaB (iaA .size());
845 ArrayRCP<LO> jaB (jaA .size());
846 ArrayRCP<SC> valB(valA.size());
847 for (
int i = 0; i < iaA .size(); i++) iaB [i] = iaA[i];
848 for (
int i = 0; i < jaA .size(); i++) jaB [i] = jaA[i];
849 for (
int i = 0; i < valA.size(); i++) valB[i] = Teuchos::ScalarTraits<SC>::magnitude(valA[i]);
851 RCP<Matrix> B = rcp(
new CrsMatrixWrap(A.getRowMap(), A.getColMap(), 0, Xpetra::StaticProfile));
852 RCP<CrsMatrix> Bcrs = rcp_dynamic_cast<CrsMatrixWrap>(B)->getCrsMatrix();
853 Bcrs->setAllValues(iaB, jaB, valB);
854 Bcrs->expertStaticFillComplete(A.getDomainMap(), A.getRangeMap());
860 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
862 return "Thyra::MueLuTpetraQ2Q1PreconditionerFactory";
868 #endif // ifdef THYRA_MUELU_TPETRA_Q2Q1PRECONDITIONER_FACTORY_DEF_HPP Generic Smoother Factory for generating the smoothers of the MG hierarchy.
void SetDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, const ParameterList ¶mList) const
This class specifies the default factory that should generate some data on a Level if the data does n...
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.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
Factory for building blocked, segregated prolongation operators.
Class that encapsulates external library smoothers.
Factory for building coarse matrices.
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void uninitializePrec(PreconditionerBase< SC > *prec, Teuchos::RCP< const LinearOpSourceBase< SC > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Base class for smoother prototypes.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Get factory associated with a particular data name.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< SC > > &fwdOp, PreconditionerBase< SC > *prec, const ESupportSolveUse supportSolveUse) const
Factory for building restriction operators using a prolongator factory.
std::string description() const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Class that holds all level-specific information.
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
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
Factory for building a thresholded operator.
AmalgamationFactory for subblocks of strided map based amalgamation data.
Factory for building the constraint operator.
direct solver for nxn blocked matrices
bool isCompatible(const LinearOpSourceBase< SC > &fwdOp) 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.
MueLu representation of a graph.
void SetBlockDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, LO row, LO col, const std::string &mode, const ParameterList ¶mList) const
Factory for building the Schur Complement for a 2x2 block matrix.
Factory for creating a graph base on a given matrix.
Teuchos::RCP< Xpetra::Matrix< SC, LO, GO, NO > > Absolute(const Xpetra::Matrix< SC, LO, GO, NO > &A) 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.
Exception throws to report errors in the internal logical of the program.
Factory for building filtered matrices using filtered graphs.
#define MUELU_GPD(name, type, defaultValue)
Teuchos::RCP< PreconditionerBase< SC > > createPrec() 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()
RCP< MueLu::FactoryBase > GetSmoother(const std::string &type, const ParameterList ¶mList, bool coarseSolver) const
static const RCP< const NoFactory > getRCP()
Static Get() functions.