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);
268 const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& velCoords,
269 const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& presCoords,
270 const ArrayRCP<LocalOrdinal>& p2vMap,
271 const Teko::LinearOp& thA11,
const Teko::LinearOp& thA12,
const Teko::LinearOp& thA21,
const Teko::LinearOp& thA11_9Pt)
const
275 typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TP_Crs;
276 typedef Tpetra::Operator <SC,LO,GO,NO> TP_Op;
278 typedef Xpetra::BlockedCrsMatrix <SC,LO,GO,NO> BlockedCrsMatrix;
279 typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
280 typedef Xpetra::CrsMatrixWrap <SC,LO,GO,NO> CrsMatrixWrap;
281 typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
282 typedef Xpetra::MapExtractor <SC,LO,GO,NO> MapExtractor;
283 typedef Xpetra::Map <LO,GO,NO> Map;
284 typedef Xpetra::MapFactory <LO,GO,NO> MapFactory;
285 typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
286 typedef Xpetra::MatrixFactory <SC,LO,GO,NO> MatrixFactory;
287 typedef Xpetra::StridedMapFactory <LO,GO,NO> StridedMapFactory;
289 typedef MueLu::Hierarchy <SC,LO,GO,NO> Hierarchy;
291 const RCP<const Teuchos::Comm<int> > comm = velCoords->getMap()->getComm();
294 RCP<Thyra::LinearOpBase<SC> > ThNonConstA11 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11);
295 RCP<Thyra::LinearOpBase<SC> > ThNonConstA21 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA21);
296 RCP<Thyra::LinearOpBase<SC> > ThNonConstA12 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA12);
297 RCP<Thyra::LinearOpBase<SC> > ThNonConstA11_9Pt = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11_9Pt);
299 RCP<TP_Op> TpetA11 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11);
300 RCP<TP_Op> TpetA21 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA21);
301 RCP<TP_Op> TpetA12 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA12);
302 RCP<TP_Op> TpetA11_9Pt = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11_9Pt);
304 RCP<TP_Crs> TpetCrsA11 = rcp_dynamic_cast<TP_Crs>(TpetA11);
305 RCP<TP_Crs> TpetCrsA21 = rcp_dynamic_cast<TP_Crs>(TpetA21);
306 RCP<TP_Crs> TpetCrsA12 = rcp_dynamic_cast<TP_Crs>(TpetA12);
307 RCP<TP_Crs> TpetCrsA11_9Pt = rcp_dynamic_cast<TP_Crs>(TpetA11_9Pt);
314 Xpetra::global_size_t numVel = A_11->getRowMap()->getLocalNumElements();
315 Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getLocalNumElements();
319 RCP<const Map> domainMap2 = tmp_A_12->getDomainMap();
320 RCP<const Map> rangeMap2 = tmp_A_21->getRangeMap();
321 Xpetra::global_size_t numRows2 = rangeMap2->getLocalNumElements();
322 Xpetra::global_size_t numCols2 = domainMap2->getLocalNumElements();
323 ArrayView<const GO> rangeElem2 = rangeMap2->getLocalElementList();
324 ArrayView<const GO> domainElem2 = domainMap2->getLocalElementList();
325 ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getLocalElementList();
326 ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getLocalElementList();
328 Xpetra::UnderlyingLib lib = domainMap2->lib();
329 GO indexBase = domainMap2->getIndexBase();
331 Array<GO> newRowElem2(numRows2, 0);
332 for (Xpetra::global_size_t i = 0; i < numRows2; i++)
333 newRowElem2[i] = numVel + rangeElem2[i];
335 RCP<const Map> newRangeMap2 = MapFactory::Build(lib, numRows2, newRowElem2, indexBase, comm);
338 Array<GO> newColElem2(numCols2, 0);
339 for (Xpetra::global_size_t i = 0; i < numCols2; i++)
340 newColElem2[i] = numVel + domainElem2[i];
342 RCP<const Map> newDomainMap2 = MapFactory::Build(lib, numCols2, newColElem2, indexBase, comm);
344 RCP<Matrix> A_12 = MatrixFactory::Build(tmp_A_12->getRangeMap(), newDomainMap2, tmp_A_12->getLocalMaxNumRowEntries());
345 RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getLocalMaxNumRowEntries());
347 RCP<CrsMatrix> A_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_11) ->getCrsMatrix();
348 RCP<CrsMatrix> A_12_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_12) ->getCrsMatrix();
349 RCP<CrsMatrix> A_21_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_21) ->getCrsMatrix();
350 RCP<CrsMatrix> A_11_crs_9Pt = rcp_dynamic_cast<CrsMatrixWrap>(A_11_9Pt)->getCrsMatrix();
353 RCP<Matrix> A_22 = MatrixFactory::Build(newRangeMap2, newDomainMap2, 1);
354 RCP<CrsMatrix> A_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_22) ->getCrsMatrix();
357 Array<SC> smallVal(1, 1.0e-10);
362 ArrayView<const LO> inds;
363 ArrayView<const SC> vals;
364 for (
LO row = 0; row < as<LO>(numRows2); ++row) {
365 tmp_A_21->getLocalRowView(row, inds, vals);
367 size_t nnz = inds.size();
368 Array<GO> newInds(nnz, 0);
369 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
370 newInds[colID] = colElem1[inds[colID]];
372 A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
373 A_22_crs->insertGlobalValues(newRowElem2[row], Array<LO>(1, newRowElem2[row]), smallVal);
375 A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
376 A_22_crs->fillComplete(newDomainMap2, newRangeMap2);
378 RCP<Matrix> A_22 = Teuchos::null;
379 RCP<CrsMatrix> A_22_crs = Teuchos::null;
381 ArrayView<const LO> inds;
382 ArrayView<const SC> vals;
383 for (
LO row = 0; row < as<LO>(numRows2); ++row) {
384 tmp_A_21->getLocalRowView(row, inds, vals);
386 size_t nnz = inds.size();
387 Array<GO> newInds(nnz, 0);
388 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
389 newInds[colID] = colElem1[inds[colID]];
391 A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
393 A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
398 for (
LO row = 0; row < as<LO>(tmp_A_12->getRowMap()->getLocalNumElements()); ++row) {
399 tmp_A_12->getLocalRowView(row, inds, vals);
401 size_t nnz = inds.size();
402 Array<GO> newInds(nnz, 0);
403 for (
LO colID = 0; colID < as<LO>(nnz); colID++)
404 newInds[colID] = newColElem2[inds[colID]];
406 A_12_crs->insertGlobalValues(rowElem1[row], newInds, vals);
408 A_12_crs->fillComplete(newDomainMap2, tmp_A_12->getRangeMap());
410 RCP<Matrix> A_12_abs = Absolute(*A_12);
411 RCP<Matrix> A_21_abs = Absolute(*A_21);
416 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
417 Teuchos::FancyOStream& out = *fancy;
418 out.setOutputToRootOnly(0);
419 RCP<Matrix> BBt = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21,
false, *A_12,
false, out);
420 RCP<Matrix> BBt_abs = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21_abs,
false, *A_12_abs,
false, out);
422 SC dropTol = (paramList.get<
int>(
"useFilters") ? paramList.get<
double>(
"tau_1") : 0.00);
423 RCP<Matrix> filteredA = FilterMatrix(*A_11, *A_11, dropTol);
424 RCP<Matrix> filteredB = FilterMatrix(*BBt, *BBt_abs, dropTol);
426 RCP<Matrix> fA_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredA);
427 RCP<Matrix> fA_12_crs = Teuchos::null;
428 RCP<Matrix> fA_21_crs = Teuchos::null;
429 RCP<Matrix> fA_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredB);
432 std::vector<size_t> stridingInfo(1, 1);
433 int stridedBlockId = -1;
435 Array<GO> elementList(numVel+numPres);
436 Array<GO> velElem = A_12_crs->getRangeMap()->getLocalElementList();
437 Array<GO> presElem = A_21_crs->getRangeMap()->getLocalElementList();
439 for (Xpetra::global_size_t i = 0 ; i < numVel; i++) elementList[i] = velElem[i];
440 for (Xpetra::global_size_t i = numVel; i < numVel+numPres; i++) elementList[i] = presElem[i-numVel];
441 RCP<const Map> fullMap = StridedMapFactory::Build(Xpetra::UseTpetra, numVel+numPres, elementList(), indexBase, stridingInfo, comm);
443 std::vector<RCP<const Map> > partMaps(2);
444 partMaps[0] = StridedMapFactory::Build(Xpetra::UseTpetra, numVel, velElem, indexBase, stridingInfo, comm);
445 partMaps[1] = StridedMapFactory::Build(Xpetra::UseTpetra, numPres, presElem, indexBase, stridingInfo, comm, stridedBlockId, numVel);
448 RCP<const MapExtractor> mapExtractor = MapExtractorFactory::Build(fullMap, partMaps);
449 RCP<BlockedCrsMatrix> fA = rcp(
new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
450 fA->setMatrix(0, 0, fA_11_crs);
451 fA->setMatrix(0, 1, fA_12_crs);
452 fA->setMatrix(1, 0, fA_21_crs);
453 fA->setMatrix(1, 1, fA_22_crs);
460 SetDependencyTree(M, paramList);
462 RCP<Hierarchy> H = rcp(
new Hierarchy);
463 RCP<MueLu::Level> finestLevel = H->GetLevel(0);
464 finestLevel->Set(
"A", rcp_dynamic_cast<Matrix>(fA));
465 finestLevel->Set(
"p2vMap", p2vMap);
466 finestLevel->Set(
"CoordinatesVelocity", Xpetra::toXpetra(velCoords));
467 finestLevel->Set(
"CoordinatesPressure", Xpetra::toXpetra(presCoords));
468 finestLevel->Set(
"AForPat", A_11_9Pt);
469 H->SetMaxCoarseSize(
MUELU_GPD(
"coarse: max size",
int, 1));
479 H->Keep(
"Ptent", M.
GetFactory(
"Ptent").get());
480 H->Setup(M, 0,
MUELU_GPD(
"max levels",
int, 3));
483 for (
int i = 1; i < H->GetNumLevels(); i++) {
484 RCP<Matrix> P = H->GetLevel(i)->template Get<RCP<Matrix> >(
"P");
485 RCP<BlockedCrsMatrix> Pcrs = rcp_dynamic_cast<BlockedCrsMatrix>(P);
486 RCP<Matrix> Pp = Pcrs->getMatrix(1,1);
487 RCP<Matrix> Pv = Pcrs->getMatrix(0,0);
489 Xpetra::IO<SC,LO,GO,NO>::Write(
"Pp_l" +
MueLu::toString(i) +
".mm", *Pp);
490 Xpetra::IO<SC,LO,GO,NO>::Write(
"Pv_l" +
MueLu::toString(i) +
".mm", *Pv);
497 std::string smootherType =
MUELU_GPD(
"smoother: type", std::string,
"vanka");
498 ParameterList smootherParams;
499 if (paramList.isSublist(
"smoother: params"))
500 smootherParams = paramList.sublist(
"smoother: params");
501 M.
SetFactory(
"Smoother", GetSmoother(smootherType, smootherParams,
false));
503 std::string coarseType =
MUELU_GPD(
"coarse: type", std::string,
"direct");
504 ParameterList coarseParams;
505 if (paramList.isSublist(
"coarse: params"))
506 coarseParams = paramList.sublist(
"coarse: params");
507 M.
SetFactory(
"CoarseSolver", GetSmoother(coarseType, coarseParams,
true));
509#ifdef HAVE_MUELU_DEBUG
513 RCP<BlockedCrsMatrix> A = rcp(
new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
514 A->setMatrix(0, 0, A_11);
515 A->setMatrix(0, 1, A_12);
516 A->setMatrix(1, 0, A_21);
517 A->setMatrix(1, 1, A_22);
520 H->GetLevel(0)->Set(
"A", rcp_dynamic_cast<Matrix>(A));
522 H->Setup(M, 0, H->GetNumLevels());
749 GetSmoother(
const std::string& type,
const ParameterList& paramList,
bool coarseSolver)
const {
750 typedef Teuchos::ParameterEntry ParameterEntry;
752 typedef MueLu::BlockedDirectSolver <SC,LO,GO,NO> BlockedDirectSolver;
753 typedef MueLu::BraessSarazinSmoother <SC,LO,GO,NO> BraessSarazinSmoother;
754 typedef MueLu::DirectSolver <SC,LO,GO,NO> DirectSolver;
755 typedef MueLu::FactoryManager <SC,LO,GO,NO> FactoryManager;
757 typedef MueLu::SmootherFactory <SC,LO,GO,NO> SmootherFactory;
758 typedef MueLu::SmootherPrototype <SC,LO,GO,NO> SmootherPrototype;
759 typedef MueLu::TrilinosSmoother <SC,LO,GO,NO> TrilinosSmoother;
761 RCP<SmootherPrototype> smootherPrototype;
762 if (type ==
"none") {
763 return Teuchos::null;
765 }
else if (type ==
"vanka") {
767 ParameterList schwarzList;
768 schwarzList.set(
"schwarz: overlap level", as<int>(0));
769 schwarzList.set(
"schwarz: zero starting solution",
false);
770 schwarzList.set(
"subdomain solver name",
"Block_Relaxation");
772 ParameterList& innerSolverList = schwarzList.sublist(
"subdomain solver parameters");
773 innerSolverList.set(
"partitioner: type",
"user");
774 innerSolverList.set(
"partitioner: overlap",
MUELU_GPD(
"partitioner: overlap",
int, 1));
775 innerSolverList.set(
"relaxation: type",
MUELU_GPD(
"relaxation: type", std::string,
"Gauss-Seidel"));
776 innerSolverList.set(
"relaxation: sweeps",
MUELU_GPD(
"relaxation: sweeps",
int, 1));
777 innerSolverList.set(
"relaxation: damping factor",
MUELU_GPD(
"relaxation: damping factor",
double, 0.5));
778 innerSolverList.set(
"relaxation: zero starting solution",
false);
781 std::string ifpackType =
"SCHWARZ";
783 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, schwarzList));
785 }
else if (type ==
"schwarz") {
787 std::string ifpackType =
"SCHWARZ";
789 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, paramList));
791 }
else if (type ==
"braess-sarazin") {
794 bool lumping =
MUELU_GPD(
"bs: lumping",
bool,
false);
796 RCP<SchurComplementFactory> schurFact = rcp(
new SchurComplementFactory());
797 schurFact->SetParameter(
"omega", ParameterEntry(omega));
798 schurFact->SetParameter(
"lumping", ParameterEntry(lumping));
802 RCP<SmootherPrototype> schurSmootherPrototype;
803 std::string schurSmootherType = (paramList.isParameter(
"schur smoother: type") ? paramList.get<std::string>(
"schur smoother: type") :
"RELAXATION");
804 if (schurSmootherType ==
"RELAXATION") {
805 ParameterList schurSmootherParams = paramList.sublist(
"schur smoother: params");
807 schurSmootherPrototype = rcp(
new TrilinosSmoother(schurSmootherType, schurSmootherParams));
809 schurSmootherPrototype = rcp(
new DirectSolver());
811 schurSmootherPrototype->SetFactory(
"A", schurFact);
813 RCP<SmootherFactory> schurSmootherFact = rcp(
new SmootherFactory(schurSmootherPrototype));
816 RCP<FactoryManager> braessManager = rcp(
new FactoryManager());
817 braessManager->SetFactory(
"A", schurFact);
818 braessManager->SetFactory(
"Smoother", schurSmootherFact);
819 braessManager->SetFactory(
"PreSmoother", schurSmootherFact);
820 braessManager->SetFactory(
"PostSmoother", schurSmootherFact);
821 braessManager->SetIgnoreUserData(
true);
823 smootherPrototype = rcp(
new BraessSarazinSmoother());
824 smootherPrototype->SetParameter(
"Sweeps", ParameterEntry(
MUELU_GPD(
"bs: sweeps",
int, 1)));
825 smootherPrototype->SetParameter(
"lumping", ParameterEntry(lumping));
826 smootherPrototype->SetParameter(
"Damping factor", ParameterEntry(omega));
827 smootherPrototype->SetParameter(
"q2q1 mode", ParameterEntry(
true));
828 rcp_dynamic_cast<BraessSarazinSmoother>(smootherPrototype)->AddFactoryManager(braessManager, 0);
830 }
else if (type ==
"ilu") {
831 std::string ifpackType =
"RILUK";
833 smootherPrototype = rcp(
new TrilinosSmoother(ifpackType, paramList));
835 }
else if (type ==
"direct") {
836 smootherPrototype = rcp(
new BlockedDirectSolver());
842 return coarseSolver ? rcp(
new SmootherFactory(smootherPrototype, Teuchos::null)) : rcp(
new SmootherFactory(smootherPrototype));