45#ifndef MUELU_MAXWELL1_DEF_HPP
46#define MUELU_MAXWELL1_DEF_HPP
52#include "Xpetra_Map.hpp"
53#include "Xpetra_CrsMatrixUtils.hpp"
54#include "Xpetra_MatrixUtils.hpp"
57#include "MueLu_Maxwell_Utils.hpp"
59#include "MueLu_ReitzingerPFactory.hpp"
60#include "MueLu_SaPFactory.hpp"
61#include "MueLu_AggregationExportFactory.hpp"
62#include "MueLu_Utilities.hpp"
64#include "MueLu_Hierarchy.hpp"
65#include "MueLu_RAPFactory.hpp"
67#include "MueLu_PerfUtils.hpp"
68#include "MueLu_ParameterListInterpreter.hpp"
70#include <MueLu_HierarchyUtils.hpp>
71# include "MueLu_Utilities_kokkos.hpp"
75#include <MueLu_RefMaxwellSmoother.hpp>
78#include "cuda_profiler_api.h"
85 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 return SM_Matrix_->getDomainMap();
91 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
93 return SM_Matrix_->getRangeMap();
97 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 if (list.isType<std::string>(
"parameterlist: syntax") && list.get<std::string>(
"parameterlist: syntax") ==
"ml") {
101 list.remove(
"parameterlist: syntax");
102 Teuchos::ParameterList newList;
109 newList.sublist(
"maxwell1: 22list").set(
"use kokkos refactor",
false);
110 newList.sublist(
"maxwell1: 22list").set(
"tentative: constant column sums",
false);
111 newList.sublist(
"maxwell1: 22list").set(
"tentative: calculate qr",
false);
113 newList.sublist(
"maxwell1: 11list").set(
"use kokkos refactor",
false);
114 newList.sublist(
"maxwell1: 11list").set(
"tentative: constant column sums",
false);
115 newList.sublist(
"maxwell1: 11list").set(
"tentative: calculate qr",
false);
117 if(list.isParameter(
"aggregation: damping factor") && list.get<
double>(
"aggregation: damping factor") == 0.0)
118 newList.sublist(
"maxwell1: 11list").set(
"multigrid algorithm",
"unsmoothed reitzinger");
120 newList.sublist(
"maxwell1: 11list").set(
"multigrid algorithm",
"smoothed reitzinger");
121 newList.sublist(
"maxwell1: 11list").set(
"aggregation: type",
"uncoupled");
123 newList.sublist(
"maxwell1: 22list").set(
"multigrid algorithm",
"unsmoothed");
124 newList.sublist(
"maxwell1: 22list").set(
"aggregation: type",
"uncoupled");
126 if (newList.sublist(
"maxwell1: 22list").isType<std::string>(
"verbosity"))
127 newList.set(
"verbosity", newList.sublist(
"maxwell1: 22list").get<std::string>(
"verbosity"));
130 std::vector<std::string> convert = {
"coarse:",
"smoother:",
"smoother: pre",
"smoother: post"};
131 for (
auto it = convert.begin(); it != convert.end(); ++it) {
132 if (newList.sublist(
"maxwell1: 22list").isType<std::string>(*it +
" type")) {
133 newList.sublist(
"maxwell1: 11list").set(*it+
" type", newList.sublist(
"maxwell1: 22list").get<std::string>(*it+
" type"));
134 newList.sublist(
"maxwell1: 22list").remove(*it+
" type");
136 if (newList.sublist(
"maxwell1: 22list").isSublist(*it+
" params")) {
137 newList.sublist(
"maxwell1: 11list").set(*it+
" params", newList.sublist(
"maxwell1: 22list").sublist(*it+
" params"));
138 newList.sublist(
"maxwell1: 22list").remove(*it+
" params");
142 newList.sublist(
"maxwell1: 22list").set(
"smoother: type",
"none");
143 newList.sublist(
"maxwell1: 22list").set(
"coarse: type",
"none");
147 std::string mode_string = list.get(
"maxwell1: mode", MasterList::getDefault<std::string>(
"maxwell1: mode"));
148 applyBCsTo22_ = list.get(
"maxwell1: apply BCs to 22",
false);
149 dump_matrices_ = list.get(
"maxwell1: dump matrices", MasterList::getDefault<bool>(
"maxwell1: dump matrices"));
152 Teuchos::ParameterList defaultSmootherList;
153 defaultSmootherList.set(
"smoother: type",
"CHEBYSHEV");
154 defaultSmootherList.sublist(
"smoother: params").set(
"chebyshev: degree",2);
155 defaultSmootherList.sublist(
"smoother: params").set(
"chebyshev: ratio eigenvalue",7.0);
156 defaultSmootherList.sublist(
"smoother: params").set(
"chebyshev: eigenvalue max iterations",30);
159 std::string verbosity = list.get(
"verbosity",
"high");
163 if(mode_ != MODE_GMHD_STANDARD) {
164 if(mode_string ==
"standard") mode_ = MODE_STANDARD;
165 else if(mode_string ==
"refmaxwell") mode_ = MODE_REFMAXWELL;
166 else if(mode_string ==
"edge only") mode_ = MODE_EDGE_ONLY;
168 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"Must use mode 'standard', 'refmaxwell', 'edge only', or use the GMHD constructor.");
174 if(list.isSublist(
"maxwell1: 22list"))
175 precList22_ = list.sublist(
"maxwell1: 22list");
176 else if(list.isSublist(
"refmaxwell: 22list"))
177 precList22_ = list.sublist(
"refmaxwell: 22list");
178 if(mode_ == MODE_EDGE_ONLY || mode_ == MODE_STANDARD || mode_ == MODE_GMHD_STANDARD)
179 precList22_.set(
"smoother: pre or post",
"none");
180 else if(!precList22_.isType<std::string>(
"Preconditioner Type") &&
181 !precList22_.isType<std::string>(
"smoother: type") &&
182 !precList22_.isType<std::string>(
"smoother: pre type") &&
183 !precList22_.isType<std::string>(
"smoother: post type")) {
184 precList22_ = defaultSmootherList;
186 precList22_.set(
"verbosity",precList22_.get(
"verbosity",verbosity));
192 if(list.isSublist(
"maxwell1: 11list"))
193 precList11_ = list.sublist(
"maxwell1: 11list");
194 else if(list.isSublist(
"refmaxwell: 11list"))
195 precList11_ = list.sublist(
"refmaxwell: 11list");
197 if(mode_ == MODE_GMHD_STANDARD) {
198 precList11_.set(
"smoother: pre or post",
"none");
199 precList11_.set(
"smoother: type",
"none");
201 if(!precList11_.isType<std::string>(
"Preconditioner Type") &&
202 !precList11_.isType<std::string>(
"smoother: type") &&
203 !precList11_.isType<std::string>(
"smoother: pre type") &&
204 !precList11_.isType<std::string>(
"smoother: post type")) {
205 if(mode_ == MODE_EDGE_ONLY || mode_ == MODE_REFMAXWELL) {
206 precList11_ = defaultSmootherList;
208 if (mode_ == MODE_STANDARD) {
209 precList11_.set(
"smoother: type",
"HIPTMAIR");
210 precList11_.sublist(
"hiptmair: smoother type 1",
"CHEBYSHEV");
211 precList11_.sublist(
"hiptmair: smoother type 2",
"CHEBYSHEV");
212 precList11_.sublist(
"hiptmair: smoother list 1") = defaultSmootherList;
213 precList11_.sublist(
"hiptmair: smoother list 2") = defaultSmootherList;
216 precList11_.set(
"verbosity",precList11_.get(
"verbosity",verbosity));
220 !precList11_.isType<std::string>(
"Preconditioner Type") &&
221 !precList11_.isParameter(
"reuse: type"))
222 precList11_.set(
"reuse: type",
"full");
224 !precList22_.isType<std::string>(
"Preconditioner Type") &&
225 !precList22_.isParameter(
"reuse: type"))
226 precList22_.set(
"reuse: type",
"full");
230# ifdef HAVE_MUELU_SERIAL
231 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosSerialWrapperNode).name())
234# ifdef HAVE_MUELU_OPENMP
235 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosOpenMPWrapperNode).name())
238# ifdef HAVE_MUELU_CUDA
239 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosCudaWrapperNode).name())
242 useKokkos_ = list.get(
"use kokkos refactor",useKokkos_);
244 parameterList_ = list;
248 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
250 Teuchos::ParameterList precListGmhd;
254 HierarchyGmhd_->GetLevel(0)->Set(
"A", GmhdA_Matrix_);
255 GmhdA_Matrix_->setObjectLabel(
"GmhdA");
257 TEUCHOS_TEST_FOR_EXCEPTION( !List.isSublist(
"maxwell1: Gmhdlist"),
Exceptions::RuntimeError,
"Must provide maxwell1: Gmhdlist for GMHD setup");
258 precListGmhd = List.sublist(
"maxwell1: Gmhdlist");
259 precListGmhd.set(
"coarse: max size",1);
260 precListGmhd.set(
"max levels",HierarchyGmhd_->GetNumLevels());
262 HierarchyGmhd_->setlib(GmhdA_Matrix_->getDomainMap()->lib());
263 HierarchyGmhd_->SetProcRankVerbose(GmhdA_Matrix_->getDomainMap()->getComm()->getRank());
264 mueLuFactory->SetupHierarchy(*HierarchyGmhd_);
267 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
282#ifdef HAVE_MUELU_CUDA
283 if (parameterList_.get<
bool>(
"maxwell1: cuda profile setup",
false)) cudaProfilerStart();
286 std::string timerLabel;
288 timerLabel =
"MueLu Maxwell1: compute (reuse)";
290 timerLabel =
"MueLu Maxwell1: compute";
291 RCP<Teuchos::TimeMonitor> tmCompute = getTimer(timerLabel);
295 bool have_generated_Kn =
false;
296 if(Kn_Matrix_.is_null()) {
297 GetOStream(
Runtime0) <<
"Maxwell1::compute(): Kn not provided. Generating." << std::endl;
298 Kn_Matrix_ = generate_kn();
299 have_generated_Kn =
true;
302 if (parameterList_.get<
bool>(
"rap: fix zero diagonals",
true)) {
304 if (parameterList_.isType<
magnitudeType>(
"rap: fix zero diagonals threshold"))
305 threshold = parameterList_.get<
magnitudeType>(
"rap: fix zero diagonals threshold",
306 Teuchos::ScalarTraits<double>::eps());
308 threshold = Teuchos::as<magnitudeType>(parameterList_.get<
double>(
"rap: fix zero diagonals threshold",
309 Teuchos::ScalarTraits<double>::eps()));
310 Scalar replacement = Teuchos::as<Scalar>(parameterList_.get<
double>(
"rap: fix zero diagonals replacement",
311 MasterList::getDefault<double>(
"rap: fix zero diagonals replacement")));
312 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Kn_Matrix_,
true, GetOStream(
Warnings1), threshold, replacement);
318 Kn_Matrix_->setObjectLabel(
"Maxwell1 (2,2)");
321 if (!Coords_.is_null())
322 precList22_.sublist(
"user data").set(
"Coordinates",Coords_);
327 if(precList22_.isParameter(
"repartition: enable")) {
328 bool repartition = precList22_.get<
bool>(
"repartition: enable");
329 precList11_.set(
"repartition: enable",repartition);
333 precList22_.set(
"repartition: rebalance P and R",
true);
335 if (precList22_.isParameter(
"repartition: use subcommunicators")) {
336 precList11_.set(
"repartition: use subcommunicators", precList22_.get<
bool>(
"repartition: use subcommunicators"));
340 if(precList11_.get<
bool>(
"repartition: use subcommunicators")==
true)
341 precList11_.set(
"repartition: use subcommunicators in place",
true);
345 precList11_.set(
"repartition: use subcommunicators",
true);
346 precList22_.set(
"repartition: use subcommunicators",
true);
350 precList11_.set(
"repartition: use subcommunicators in place",
true);
355 precList11_.remove(
"repartition: enable",
false);
376 magnitudeType rowSumTol = precList11_.get(
"aggregation: row sum drop tol",-1.0);
378 useKokkos_,BCrowsKokkos_,BCcolsKokkos_,BCdomainKokkos_,
379 BCedges_,BCnodes_,BCrows_,BCcols_,BCdomain_,
380 allEdgesBoundary_,allNodesBoundary_);
382 GetOStream(
Statistics2) <<
"MueLu::Maxwell1::compute(): Detected " << BCedges_ <<
" BC rows and " << BCnodes_ <<
" BC columns." << std::endl;
386 if (allEdgesBoundary_) {
389 GetOStream(
Warnings0) <<
"All edges are detected as boundary edges!" << std::endl;
390 mode_ = MODE_EDGE_ONLY;
393 precList22_.set(
"max levels", 1);
397 if (allNodesBoundary_) {
400 GetOStream(
Warnings0) <<
"All nodes are detected as boundary nodes!" << std::endl;
401 mode_ = MODE_EDGE_ONLY;
404 precList22_.set(
"max levels", 1);
415 D0_Matrix_->resumeFill();
417 if (D0_Matrix_->getRowMap()->lib() == Xpetra::UseEpetra)
418 replaceWith= Teuchos::ScalarTraits<SC>::eps();
420 replaceWith = Teuchos::ScalarTraits<SC>::zero();
423 GetOStream(
Runtime0) <<
"Maxwell1::compute(): nuking BC rows/cols of D0" << std::endl;
431 GetOStream(
Runtime0) <<
"Maxwell1::compute(): nuking BC rows of D0" << std::endl;
439 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
450 RCP<Matrix> Kn_Smoother_0;
451 if(have_generated_Kn) {
452 Kn_Smoother_0 = Kn_Matrix_;
455 Kn_Smoother_0 = generate_kn();
460 Hierarchy11_ = rcp(
new Hierarchy(
"Maxwell1 (1,1)"));
461 RCP<Matrix> OldSmootherMatrix;
462 for(
int i=0; i<Hierarchy22_->GetNumLevels(); i++) {
463 Hierarchy11_->AddNewLevel();
464 RCP<Level> NodeL = Hierarchy22_->GetLevel(i);
465 RCP<Level> EdgeL = Hierarchy11_->GetLevel(i);
466 RCP<Operator> NodeOp = NodeL->Get<RCP<Operator> >(
"A");
467 RCP<Matrix> NodeAggMatrix = rcp_dynamic_cast<Matrix>(NodeOp);
471 EdgeL->Set(
"A", SM_Matrix_);
472 EdgeL->Set(
"D0", D0_Matrix_);
474 EdgeL->Set(
"NodeAggMatrix",NodeAggMatrix);
475 EdgeL->Set(
"NodeMatrix",Kn_Smoother_0);
476 OldSmootherMatrix = Kn_Smoother_0;
480 auto NodalP = NodeL->Get<RCP<Matrix> >(
"P");
481 EdgeL->Set(
"Pnodal",NodalP);
485 if(!NodeAggMatrix.is_null()) {
486 EdgeL->Set(
"NodeAggMatrix",NodeAggMatrix);
487 if(!have_generated_Kn) {
490 EdgeL->Set(
"NodeMatrix",NewKn);
491 OldSmootherMatrix = NewKn;
495 EdgeL->Set(
"NodeMatrix",NodeAggMatrix);
500 EdgeL->Set(
"NodeMatrix",NodeOp);
501 EdgeL->Set(
"NodeAggMatrix",NodeOp);
507 if(NodeL->IsAvailable(
"Importer")) {
508 auto importer = NodeL->Get<RCP<const Import> >(
"Importer");
509 EdgeL->Set(
"NodeImporter",importer);
516 std::string fine_smoother = precList11_.get<std::string>(
"smoother: type");
518 SM_Matrix_->setObjectLabel(
"A(1,1)");
519 precList11_.set(
"coarse: max size",1);
520 precList11_.set(
"max levels",Hierarchy22_->GetNumLevels());
522 const bool refmaxwellCoarseSolve = (precList11_.get<std::string>(
"coarse: type",
523 MasterList::getDefault<std::string>(
"coarse: type")) ==
"RefMaxwell");
524 if (refmaxwellCoarseSolve) {
525 GetOStream(
Runtime0) <<
"Maxwell1::compute(): Will set up RefMaxwell coarse solver" << std::endl;
526 precList11_.set(
"coarse: type",
"none");
530 Teuchos::ParameterList nonSerialList11, processedPrecList11;
533 Hierarchy11_->setlib(SM_Matrix_->getDomainMap()->lib());
534 Hierarchy11_->SetProcRankVerbose(SM_Matrix_->getDomainMap()->getComm()->getRank());
537 mueLuFactory->SetupHierarchy(*Hierarchy11_);
539 if (refmaxwellCoarseSolve) {
540 GetOStream(
Runtime0) <<
"Maxwell1::compute(): Setting up RefMaxwell coarse solver" << std::endl;
541 auto coarseLvl = Hierarchy11_->GetLevel(Hierarchy11_->GetNumLevels()-1);
543 precList11_.sublist(
"coarse: params")));
544 coarseSolver->Setup(*coarseLvl);
545 coarseLvl->Set(
"PreSmoother",
546 rcp_dynamic_cast<SmootherBase>(coarseSolver,
true));
549 if(mode_ == MODE_REFMAXWELL) {
550 if(Hierarchy11_->GetNumLevels() > 1) {
551 RCP<Level> EdgeL = Hierarchy11_->GetLevel(1);
552 P11_ = EdgeL->Get<RCP<Matrix> >(
"P");
564#ifdef MUELU_MAXWELL1_DEBUG
565 for(
int i=0; i<Hierarchy11_->GetNumLevels(); i++) {
566 RCP<Level> L = Hierarchy11_->GetLevel(i);
567 RCP<Matrix> EdgeMatrix = rcp_dynamic_cast<Matrix>(L->Get<RCP<Operator> >(
"A"));
568 RCP<Matrix> NodeMatrix = rcp_dynamic_cast<Matrix>(L->Get<RCP<Operator> >(
"NodeMatrix"));
569 RCP<Matrix> NodeAggMatrix = rcp_dynamic_cast<Matrix>(L->Get<RCP<Operator> >(
"NodeAggMatrix"));
570 RCP<Matrix> D0 =rcp_dynamic_cast<Matrix>( L->Get<RCP<Operator> >(
"D0"));
572 auto nrmE = EdgeMatrix->getFrobeniusNorm();
573 auto nrmN = NodeMatrix->getFrobeniusNorm();
574 auto nrmNa = NodeAggMatrix->getFrobeniusNorm();
575 auto nrmD0= D0->getFrobeniusNorm();
577 std::cout<<
"DEBUG: Norms on Level "<<i<<
" E/N/NA/D0 = "<<nrmE<<
" / "<<nrmN <<
" / "<<nrmNa<<
" / "<< nrmD0 <<std::endl;
578 std::cout<<
"DEBUG: NNZ on Level "<<i<<
" E/N/NA/D0 = "<<
579 EdgeMatrix->getGlobalNumEntries()<<
" / "<<
580 NodeMatrix->getGlobalNumEntries()<<
" / "<<
581 NodeAggMatrix->getGlobalNumEntries()<<
" / "<<
582 D0->getGlobalNumEntries()<<std::endl;
589 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
595 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu Maxwell1: Build Kn");
597 Level fineLevel, coarseLevel;
603 fineLevel.
Set(
"A",SM_Matrix_);
604 coarseLevel.
Set(
"P",D0_Matrix_);
607 coarseLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
608 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
609 coarseLevel.setObjectLabel(
"Maxwell1 (2,2)");
610 fineLevel.setObjectLabel(
"Maxwell1 (2,2)");
612 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
613 ParameterList rapList = *(rapFact->GetValidParameterList());
614 rapList.set(
"transpose: use implicit",
true);
615 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
616 rapFact->SetParameterList(rapList);
617 coarseLevel.
Request(
"A", rapFact.get());
619 coarseLevel.
Request(
"AP reuse data", rapFact.get());
620 coarseLevel.
Request(
"RAP reuse data", rapFact.get());
623 RCP<Matrix> Kn_Matrix = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
624 Kn_Matrix->setObjectLabel(
"A(2,2)");
631 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
633 if(mode_ == MODE_REFMAXWELL) {
634 RCP<Teuchos::TimeMonitor> tmAlloc = getTimer(
"MueLu Maxwell1: Allocate MVs");
636 residualFine_ = MultiVectorFactory::Build(SM_Matrix_->getRangeMap(), numVectors);
638 if(!Hierarchy11_.is_null() && Hierarchy11_->GetNumLevels() > 1) {
639 RCP<Level> EdgeL = Hierarchy11_->GetLevel(1);
640 RCP<Matrix> A = EdgeL->Get<RCP<Matrix> >(
"A");
641 residual11c_ = MultiVectorFactory::Build(A->getRangeMap(), numVectors);
642 update11c_ = MultiVectorFactory::Build(A->getDomainMap(), numVectors);
645 if(!Hierarchy22_.is_null()) {
646 residual22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
647 update22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
655 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
657 if (dump_matrices_) {
658 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
659 Xpetra::IO<SC, LO, GO, NO>::Write(name, A);
664 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
666 if (dump_matrices_) {
667 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
668 Xpetra::IO<SC, LO, GO, NO>::Write(name, X);
673 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
675 if (dump_matrices_) {
676 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
677 Xpetra::IO<coordinateType, LO, GO, NO>::Write(name, X);
682 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
684 if (dump_matrices_) {
685 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
686 std::ofstream out(name);
687 for (
size_t i = 0; i < Teuchos::as<size_t>(v.size()); i++)
692 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
694 if (dump_matrices_) {
695 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
696 std::ofstream out(name);
697 auto vH = Kokkos::create_mirror_view (v);
698 Kokkos::deep_copy(vH , v);
699 for (
size_t i = 0; i < vH.size(); i++)
700 out << vH[i] <<
"\n";
704 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
708 return Teuchos::rcp(
new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name)));
711 return Teuchos::rcp(
new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), SM_Matrix_->getRowMap()->getComm().ptr()));
713 return Teuchos::rcp(
new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), comm.ptr()));
716 return Teuchos::null;
722 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
724 bool reuse = !SM_Matrix_.is_null();
725 SM_Matrix_ = SM_Matrix_new;
726 dump(*SM_Matrix_,
"SM.m");
732 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
735 const SC one = Teuchos::ScalarTraits<SC>::one();
736 if (!allEdgesBoundary_ && X.getNumVectors() != residualFine_->getNumVectors())
737 allocateMemory(X.getNumVectors());
739 TEUCHOS_TEST_FOR_EXCEPTION(Hierarchy11_.is_null() || Hierarchy11_->GetNumLevels() == 0,
Exceptions::RuntimeError,
"(1,1) Hiearchy is null.");
742 RCP<Level> Fine = Hierarchy11_->GetLevel(0);
743 if (Fine->IsAvailable(
"PreSmoother")) {
744 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: PreSmoother");
745 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
746 preSmoo->Apply(X, RHS,
true);
751 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: residual calculation");
756 if(!P11_.is_null()) {
757 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: (1,1) correction");
758 P11_->apply(*residualFine_,*residual11c_,Teuchos::TRANS);
759 Hierarchy11_->Iterate(*residual11c_,*update11c_,
true,1);
763 if (!allNodesBoundary_) {
764 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: (2,2) correction");
765 D0_Matrix_->apply(*residualFine_,*residual22_,Teuchos::TRANS);
766 Hierarchy22_->Iterate(*residual22_,*update22_,
true,0);
771 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: Orolongation");
773 P11_->apply(*update11c_,X,Teuchos::NO_TRANS,one,one);
774 if (!allNodesBoundary_)
775 D0_Matrix_->apply(*update22_,X,Teuchos::NO_TRANS,one,one);
779 if (Fine->IsAvailable(
"PostSmoother")) {
780 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: PostSmoother");
781 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
782 postSmoo->Apply(X, RHS,
false);
788 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
790 Hierarchy11_->Iterate(RHS,X,1,
true);
793 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
798 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu Maxwell1: solve");
799 if(mode_ == MODE_GMHD_STANDARD)
800 HierarchyGmhd_->Iterate(RHS,X,1,
true);
801 else if(mode_ == MODE_STANDARD || mode_ == MODE_EDGE_ONLY)
802 applyInverseStandard(RHS,X);
803 else if(mode_ == MODE_REFMAXWELL)
804 applyInverseRefMaxwellAdditive(RHS,X);
806 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"Must use mode 'standard', 'refmaxwell' or 'edge only' when not doing GMHD.");
811 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
817 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
819 initialize(
const Teuchos::RCP<Matrix> & D0_Matrix,
820 const Teuchos::RCP<Matrix> & Kn_Matrix,
821 const Teuchos::RCP<MultiVector> & Nullspace,
822 const Teuchos::RCP<RealValuedMultiVector> & Coords,
823 Teuchos::ParameterList& List)
826 TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
828#ifdef HAVE_MUELU_DEBUG
829 if(!Kn_Matrix.is_null()) {
830 TEUCHOS_ASSERT(Kn_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getDomainMap()));
831 TEUCHOS_ASSERT(Kn_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getDomainMap()));
835 TEUCHOS_ASSERT(D0_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getRowMap()));
838 Hierarchy11_ = Teuchos::null;
839 Hierarchy22_ = Teuchos::null;
840 HierarchyGmhd_ = Teuchos::null;
841 if (mode_ != MODE_GMHD_STANDARD) mode_ = MODE_STANDARD;
845 allEdgesBoundary_=
false;
846 allNodesBoundary_=
false;
847 dump_matrices_ =
false;
850 applyBCsTo22_ =
false;
856 if (D0_Matrix->getRowMap()->lib() == Xpetra::UseTpetra) {
861 RCP<Matrix> D0copy = MatrixFactory::Build(D0_Matrix->getRowMap(), D0_Matrix->getColMap(), 0);
862 RCP<CrsMatrix> D0copyCrs = rcp_dynamic_cast<CrsMatrixWrap>(D0copy,
true)->getCrsMatrix();
863 ArrayRCP<const size_t> D0rowptr_RCP;
864 ArrayRCP<const LO> D0colind_RCP;
865 ArrayRCP<const SC> D0vals_RCP;
866 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getAllValues(D0rowptr_RCP,
870 ArrayRCP<size_t> D0copyrowptr_RCP;
871 ArrayRCP<LO> D0copycolind_RCP;
872 ArrayRCP<SC> D0copyvals_RCP;
873 D0copyCrs->allocateAllValues(D0vals_RCP.size(),D0copyrowptr_RCP,D0copycolind_RCP,D0copyvals_RCP);
874 D0copyrowptr_RCP.deepCopy(D0rowptr_RCP());
875 D0copycolind_RCP.deepCopy(D0colind_RCP());
876 D0copyvals_RCP.deepCopy(D0vals_RCP());
877 D0copyCrs->setAllValues(D0copyrowptr_RCP,
880 D0copyCrs->expertStaticFillComplete(D0_Matrix->getDomainMap(), D0_Matrix->getRangeMap(),
881 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getImporter(),
882 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getExporter());
885 D0_Matrix_ = MatrixFactory::BuildCopy(D0_Matrix);
888 Kn_Matrix_ = Kn_Matrix;
890 Nullspace_ = Nullspace;
892 if(!Kn_Matrix_.is_null()) dump(*Kn_Matrix_,
"Kn.m");
893 if (!Nullspace_.is_null()) dump(*Nullspace_,
"nullspace.m");
894 if (!Coords_.is_null()) dumpCoords(*Coords_,
"coords.m");
900 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
902 describe(Teuchos::FancyOStream& out,
const Teuchos::EVerbosityLevel )
const {
904 std::ostringstream oss;
906 RCP<const Teuchos::Comm<int> > comm = SM_Matrix_->getDomainMap()->getComm();
910 if (!Kn_Matrix_.is_null())
911 root = comm->getRank();
916 reduceAll(*comm, Teuchos::REDUCE_MAX, root, Teuchos::ptr(&actualRoot));
921 oss <<
"\n--------------------------------------------------------------------------------\n" <<
922 "--- Maxwell1 Summary ---\n"
923 "--------------------------------------------------------------------------------" << std::endl;
929 SM_Matrix_->getRowMap()->getComm()->barrier();
931 numRows = SM_Matrix_->getGlobalNumRows();
932 nnz = SM_Matrix_->getGlobalNumEntries();
934 Xpetra::global_size_t tt = numRows;
935 int rowspacer = 3;
while (tt != 0) { tt /= 10; rowspacer++; }
937 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
939 oss <<
"block " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " << std::setw(9) <<
" nnz/row" << std::endl;
940 oss <<
"(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
942 if (!Kn_Matrix_.is_null()) {
944 numRows = Kn_Matrix_->getGlobalNumRows();
945 nnz = Kn_Matrix_->getGlobalNumEntries();
947 oss <<
"(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
953 std::string outstr = oss.str();
956 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
957 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
959 int strLength = outstr.size();
960 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
961 if (comm->getRank() != root)
962 outstr.resize(strLength);
963 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
968 if (!Hierarchy11_.is_null())
969 Hierarchy11_->describe(out, GetVerbLevel());
971 if (!Hierarchy22_.is_null())
972 Hierarchy22_->describe(out, GetVerbLevel());
974 if (!HierarchyGmhd_.is_null())
975 HierarchyGmhd_->describe(out, GetVerbLevel());
979 std::ostringstream oss2;
981 oss2 <<
"Sub-solver distribution over ranks" << std::endl;
982 oss2 <<
"( (1,1) block only is indicated by '1', (2,2) block only by '2', and both blocks by 'B' and none by '.')" << std::endl;
984 int numProcs = comm->getSize();
986 RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
988 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
992 if (!Kn_Matrix_.is_null())
994 std::vector<char> states(numProcs, 0);
996 MPI_Gather(&status, 1, MPI_CHAR, &states[0], 1, MPI_CHAR, 0, *rawMpiComm);
998 states.push_back(status);
1001 int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
1002 for (
int proc = 0; proc < numProcs; proc += rowWidth) {
1003 for (
int j = 0; j < rowWidth; j++)
1004 if (proc + j < numProcs)
1005 if (states[proc+j] == 0)
1007 else if (states[proc+j] == 1)
1009 else if (states[proc+j] == 2)
1016 oss2 <<
" " << proc <<
":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
1029#define MUELU_MAXWELL1_SHORT
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Exception throws to report errors in the internal logical of the program.
static void AddNonSerializableDataToHierarchy(HierarchyManager &HM, Hierarchy &H, const ParameterList &nonSerialList)
Add non-serializable data to Hierarchy.
static void CopyBetweenHierarchies(Hierarchy &fromHierarchy, Hierarchy &toHierarchy, const std::string fromLabel, const std::string toLabel, const std::string dataType)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Class that holds all level-specific information.
void setlib(Xpetra::UnderlyingLib lib2)
void SetLevelID(int levelID)
Set level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
void SetPreviousLevel(const RCP< Level > &previousLevel)
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
static std::string translate(Teuchos::ParameterList ¶mList, const std::string &defaultVals="")
: Translate ML parameters to MueLu parameter XML string
Xpetra::MultiVector< coordinateType, LO, GO, NO > RealValuedMultiVector
Teuchos::RCP< Matrix > generate_kn() const
Generates the Kn matrix.
Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
void compute(bool reuse=false)
Setup the preconditioner.
void GMHDSetupHierarchy(Teuchos::ParameterList &List) const
Sets up hiearchy for GMHD matrices that include generalized Ohms law equations.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
void applyInverseRefMaxwellAdditive(const MultiVector &RHS, MultiVector &X) const
apply RefMaxwell additive 2x2 style cycle
void applyInverseStandard(const MultiVector &RHS, MultiVector &X) const
apply standard Maxwell1 cycle
void initialize(const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &Kn_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List)
void dumpCoords(const RealValuedMultiVector &X, std::string name) const
dump out real-valued multivector
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
void dump(const Matrix &A, std::string name) const
dump out matrix
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
Teuchos::RCP< Teuchos::TimeMonitor > getTimer(std::string name, RCP< const Teuchos::Comm< int > > comm=Teuchos::null) const
get a (synced) timer
void setParameters(Teuchos::ParameterList &list)
Set parameters.
void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
void allocateMemory(int numVectors) const
allocate multivectors for solve
static void detectBoundaryConditionsSM(RCP< Matrix > &SM_Matrix, RCP< Matrix > &D0_Matrix, magnitudeType rowSumTol, bool useKokkos_, Kokkos::View< bool *, typename Node::device_type > &BCrowsKokkos, Kokkos::View< bool *, typename Node::device_type > &BCcolsKokkos, Kokkos::View< bool *, typename Node::device_type > &BCdomainKokkos, int &BCedges, int &BCnodes, Teuchos::ArrayRCP< bool > &BCrows, Teuchos::ArrayRCP< bool > &BCcols, Teuchos::ArrayRCP< bool > &BCdomain, bool &allEdgesBoundary, bool &allNodesBoundary)
Detect Dirichlet boundary conditions.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > PtAPWrapper(RCP< Matrix > &A, RCP< Matrix > &P, Teuchos::ParameterList ¶ms, std::string &label)
Factory for building coarse matrices.
Class that encapsulates Operator smoothers.
static void ZeroDirichletRows(RCP< Matrix > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static void ZeroDirichletCols(RCP< Matrix > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletCols, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static void ZeroDirichletCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
Namespace for MueLu classes and methods.
long ExtractNonSerializableData(const Teuchos::ParameterList &inList, Teuchos::ParameterList &serialList, Teuchos::ParameterList &nonSerialList)
Extract non-serializable data from level-specific sublists and move it to a separate parameter list.
@ Warnings0
Important warning messages (one line)
@ Statistics2
Print even more statistics.
@ Runtime0
One-liner description of what is happening.
@ Warnings1
Additional warnings.
@ Timings
Print all timing information.
MsgType toVerbLevel(const std::string &verbLevelStr)
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix,...