45 #ifndef MUELU_MAXWELL1_DEF_HPP
46 #define MUELU_MAXWELL1_DEF_HPP
53 #include "Xpetra_Map.hpp"
58 #include "MueLu_Maxwell_Utils.hpp"
60 #include "MueLu_ReitzingerPFactory.hpp"
61 #include "MueLu_Utilities.hpp"
63 #include "MueLu_Hierarchy.hpp"
64 #include "MueLu_RAPFactory.hpp"
66 #include "MueLu_PerfUtils.hpp"
67 #include "MueLu_ParameterListInterpreter.hpp"
69 #include <MueLu_HierarchyUtils.hpp>
73 #include <MueLu_RefMaxwellSmoother.hpp>
75 #ifdef HAVE_MUELU_CUDA
76 #include "cuda_profiler_api.h"
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83 return SM_Matrix_->getDomainMap();
86 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
88 return SM_Matrix_->getRangeMap();
91 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
93 if (list.
isType<std::string>(
"parameterlist: syntax") && list.
get<std::string>(
"parameterlist: syntax") ==
"ml") {
94 list.
remove(
"parameterlist: syntax");
101 newList.
sublist(
"maxwell1: 22list").
set(
"use kokkos refactor",
false);
103 newList.
sublist(
"maxwell1: 11list").
set(
"use kokkos refactor",
false);
104 newList.
sublist(
"maxwell1: 11list").
set(
"tentative: constant column sums",
false);
105 newList.
sublist(
"maxwell1: 11list").
set(
"tentative: calculate qr",
false);
107 newList.
sublist(
"maxwell1: 11list").
set(
"aggregation: use ml scaling of drop tol",
true);
108 newList.
sublist(
"maxwell1: 22list").
set(
"aggregation: use ml scaling of drop tol",
true);
110 newList.
sublist(
"maxwell1: 22list").
set(
"aggregation: min agg size", 3);
111 newList.
sublist(
"maxwell1: 22list").
set(
"aggregation: match ML phase1",
true);
112 newList.
sublist(
"maxwell1: 22list").
set(
"aggregation: match ML phase2a",
true);
113 newList.
sublist(
"maxwell1: 22list").
set(
"aggregation: match ML phase2b",
true);
115 if (list.
isParameter(
"aggregation: damping factor") && list.
get<
double>(
"aggregation: damping factor") == 0.0)
116 newList.
sublist(
"maxwell1: 11list").
set(
"multigrid algorithm",
"unsmoothed reitzinger");
118 newList.
sublist(
"maxwell1: 11list").
set(
"multigrid algorithm",
"smoothed reitzinger");
119 newList.
sublist(
"maxwell1: 11list").
set(
"aggregation: type",
"uncoupled");
121 newList.
sublist(
"maxwell1: 22list").
set(
"multigrid algorithm",
"unsmoothed");
122 newList.
sublist(
"maxwell1: 22list").
set(
"aggregation: type",
"uncoupled");
124 if (newList.
sublist(
"maxwell1: 22list").
isType<std::string>(
"verbosity"))
125 newList.
set(
"verbosity", newList.
sublist(
"maxwell1: 22list").
get<std::string>(
"verbosity"));
128 std::vector<std::string> convert = {
"coarse:",
"smoother:",
"smoother: pre",
"smoother: post"};
129 for (
auto it = convert.begin(); it != convert.end(); ++it) {
130 if (newList.
sublist(
"maxwell1: 22list").
isType<std::string>(*it +
" type")) {
131 newList.
sublist(
"maxwell1: 11list").
set(*it +
" type", newList.
sublist(
"maxwell1: 22list").
get<std::string>(*it +
" type"));
135 newList.
sublist(
"maxwell1: 11list").
set(*it +
" params", newList.
sublist(
"maxwell1: 22list").
sublist(*it +
" params"));
140 newList.
sublist(
"maxwell1: 22list").
set(
"smoother: type",
"none");
141 newList.
sublist(
"maxwell1: 22list").
set(
"coarse: type",
"none");
143 newList.
set(
"maxwell1: nodal smoother fix zero diagonal threshold", 1e-10);
144 newList.
sublist(
"maxwell1: 22list").
set(
"rap: fix zero diagonals",
true);
145 newList.
sublist(
"maxwell1: 22list").
set(
"rap: fix zero diagonals threshold", 1e-10);
149 std::string mode_string = list.
get(
"maxwell1: mode", MasterList::getDefault<std::string>(
"maxwell1: mode"));
150 applyBCsTo22_ = list.
get(
"maxwell1: apply BCs to 22",
false);
151 dump_matrices_ = list.
get(
"maxwell1: dump matrices", MasterList::getDefault<bool>(
"maxwell1: dump matrices"));
155 defaultSmootherList.
set(
"smoother: type",
"CHEBYSHEV");
156 defaultSmootherList.
sublist(
"smoother: params").
set(
"chebyshev: degree", 2);
157 defaultSmootherList.
sublist(
"smoother: params").
set(
"chebyshev: ratio eigenvalue", 7.0);
158 defaultSmootherList.
sublist(
"smoother: params").
set(
"chebyshev: eigenvalue max iterations", 30);
161 std::string verbosity = list.
get(
"verbosity",
"high");
165 if (mode_ != MODE_GMHD_STANDARD) {
166 if (mode_string ==
"standard")
167 mode_ = MODE_STANDARD;
168 else if (mode_string ==
"refmaxwell")
169 mode_ = MODE_REFMAXWELL;
170 else if (mode_string ==
"edge only")
171 mode_ = MODE_EDGE_ONLY;
180 precList22_ = list.
sublist(
"maxwell1: 22list");
181 else if (list.
isSublist(
"refmaxwell: 22list"))
182 precList22_ = list.
sublist(
"refmaxwell: 22list");
183 if (mode_ == MODE_EDGE_ONLY || mode_ == MODE_STANDARD || mode_ == MODE_GMHD_STANDARD)
184 precList22_.
set(
"smoother: pre or post",
"none");
185 else if (!precList22_.isType<std::string>(
"Preconditioner Type") &&
186 !precList22_.isType<std::string>(
"smoother: type") &&
187 !precList22_.isType<std::string>(
"smoother: pre type") &&
188 !precList22_.isType<std::string>(
"smoother: post type")) {
189 precList22_ = defaultSmootherList;
191 precList22_.
set(
"verbosity", precList22_.get(
"verbosity", verbosity));
196 precList11_ = list.
sublist(
"maxwell1: 11list");
197 else if (list.
isSublist(
"refmaxwell: 11list"))
198 precList11_ = list.
sublist(
"refmaxwell: 11list");
200 if (mode_ == MODE_GMHD_STANDARD) {
201 precList11_.
set(
"smoother: pre or post",
"none");
202 precList11_.
set(
"smoother: type",
"none");
204 if (!precList11_.isType<std::string>(
"Preconditioner Type") &&
205 !precList11_.isType<std::string>(
"smoother: type") &&
206 !precList11_.isType<std::string>(
"smoother: pre type") &&
207 !precList11_.isType<std::string>(
"smoother: post type")) {
208 if (mode_ == MODE_EDGE_ONLY || mode_ == MODE_REFMAXWELL) {
209 precList11_ = defaultSmootherList;
211 if (mode_ == MODE_STANDARD) {
212 precList11_.
set(
"smoother: type",
"HIPTMAIR");
213 precList11_.
set(
"hiptmair: smoother type 1",
"CHEBYSHEV");
214 precList11_.
set(
"hiptmair: smoother type 2",
"CHEBYSHEV");
215 precList11_.
sublist(
"hiptmair: smoother list 1") = defaultSmootherList;
216 precList11_.
sublist(
"hiptmair: smoother list 2") = defaultSmootherList;
219 precList11_.
set(
"verbosity", precList11_.get(
"verbosity", verbosity));
223 !precList11_.isType<std::string>(
"Preconditioner Type") &&
224 !precList11_.isParameter(
"reuse: type"))
225 precList11_.
set(
"reuse: type",
"full");
227 !precList22_.isType<std::string>(
"Preconditioner Type") &&
228 !precList22_.isParameter(
"reuse: type"))
229 precList22_.
set(
"reuse: type",
"full");
232 useKokkos_ = !Node::is_serial;
233 useKokkos_ = list.
get(
"use kokkos refactor", useKokkos_);
235 parameterList_ = list;
238 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
244 HierarchyGmhd_->GetLevel(0)->Set(
"A", GmhdA_Matrix_);
245 GmhdA_Matrix_->setObjectLabel(
"GmhdA");
248 precListGmhd = List.
sublist(
"maxwell1: Gmhdlist");
249 precListGmhd.
set(
"coarse: max size", 1);
250 precListGmhd.
set(
"max levels", HierarchyGmhd_->GetNumLevels());
252 HierarchyGmhd_->setlib(GmhdA_Matrix_->getDomainMap()->lib());
253 HierarchyGmhd_->SetProcRankVerbose(GmhdA_Matrix_->getDomainMap()->getComm()->getRank());
254 mueLuFactory->SetupHierarchy(*HierarchyGmhd_);
257 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
271 #ifdef HAVE_MUELU_CUDA
272 if (parameterList_.get<
bool>(
"maxwell1: cuda profile setup",
false)) cudaProfilerStart();
275 std::string timerLabel;
277 timerLabel =
"MueLu Maxwell1: compute (reuse)";
279 timerLabel =
"MueLu Maxwell1: compute";
284 bool have_generated_Kn =
false;
285 if (Kn_Matrix_.is_null()) {
287 GetOStream(
Runtime0) <<
"Maxwell1::compute(): Kn not provided. Generating." << std::endl;
288 Kn_Matrix_ = generate_kn();
289 have_generated_Kn =
true;
290 }
else if (parameterList_.get<
bool>(
"rap: fix zero diagonals",
true)) {
292 if (parameterList_.isType<
magnitudeType>(
"rap: fix zero diagonals threshold"))
293 threshold = parameterList_.get<
magnitudeType>(
"rap: fix zero diagonals threshold",
296 threshold = Teuchos::as<magnitudeType>(parameterList_.get<
double>(
"rap: fix zero diagonals threshold",
298 Scalar replacement = Teuchos::as<Scalar>(parameterList_.get<
double>(
"rap: fix zero diagonals replacement",
299 MasterList::getDefault<double>(
"rap: fix zero diagonals replacement")));
305 Kn_Matrix_->setObjectLabel(
"Maxwell1 (2,2)");
308 if (!Coords_.is_null())
309 precList22_.sublist(
"user data").set(
"Coordinates", Coords_);
314 if (precList22_.isParameter(
"repartition: enable")) {
315 bool repartition = precList22_.get<
bool>(
"repartition: enable");
316 precList11_.set(
"repartition: enable", repartition);
320 precList22_.set(
"repartition: rebalance P and R",
true);
322 if (precList22_.isParameter(
"repartition: use subcommunicators")) {
323 precList11_.set(
"repartition: use subcommunicators", precList22_.get<
bool>(
"repartition: use subcommunicators"));
327 if (precList11_.get<
bool>(
"repartition: use subcommunicators") ==
true)
328 precList11_.set(
"repartition: use subcommunicators in place",
true);
331 precList11_.set(
"repartition: use subcommunicators",
true);
332 precList22_.set(
"repartition: use subcommunicators",
true);
336 precList11_.set(
"repartition: use subcommunicators in place",
true);
340 precList11_.remove(
"repartition: enable",
false);
359 magnitudeType rowSumTol = precList11_.get(
"aggregation: row sum drop tol", -1.0);
361 useKokkos_, BCrowsKokkos_, BCcolsKokkos_, BCdomainKokkos_,
362 BCedges_, BCnodes_, BCrows_, BCcols_, BCdomain_,
363 allEdgesBoundary_, allNodesBoundary_);
365 GetOStream(
Statistics2) <<
"MueLu::Maxwell1::compute(): Detected " << BCedges_ <<
" BC rows and " << BCnodes_ <<
" BC columns." << std::endl;
369 if (allEdgesBoundary_) {
372 GetOStream(
Warnings0) <<
"All edges are detected as boundary edges!" << std::endl;
373 mode_ = MODE_EDGE_ONLY;
376 precList22_.set(
"max levels", 1);
379 if (allNodesBoundary_) {
382 GetOStream(
Warnings0) <<
"All nodes are detected as boundary nodes!" << std::endl;
383 mode_ = MODE_EDGE_ONLY;
386 precList22_.set(
"max levels", 1);
396 D0_Matrix_->resumeFill();
404 GetOStream(
Runtime0) <<
"Maxwell1::compute(): nuking BC rows/cols of D0" << std::endl;
411 GetOStream(
Runtime0) <<
"Maxwell1::compute(): nuking BC rows of D0" << std::endl;
419 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(), D0_Matrix_->getRangeMap());
430 if (have_generated_Kn) {
431 Kn_Smoother_0 = Kn_Matrix_;
433 Kn_Smoother_0 = generate_kn();
441 for (
int i = 0; i < Hierarchy22_->GetNumLevels(); i++) {
442 Hierarchy11_->AddNewLevel();
446 RCP<Matrix> NodeAggMatrix = rcp_dynamic_cast<Matrix>(NodeOp);
450 EdgeL->
Set(
"A", SM_Matrix_);
451 EdgeL->
Set(
"D0", D0_Matrix_);
453 EdgeL->
Set(
"NodeAggMatrix", NodeAggMatrix);
454 EdgeL->
Set(
"NodeMatrix", Kn_Smoother_0);
455 OldSmootherMatrix = Kn_Smoother_0;
456 OldEdgeLevel = EdgeL;
465 EdgeL->
Set(
"Pnodal", NodalP_ones);
469 if (!NodeAggMatrix.
is_null()) {
470 EdgeL->
Set(
"NodeAggMatrix", NodeAggMatrix);
471 if (!have_generated_Kn) {
482 RAPlist.
set(
"rap: fix zero diagonals",
false);
484 EdgeL->
Set(
"NodeMatrix", NewKn);
487 double thresh = parameterList_.get(
"maxwell1: nodal smoother fix zero diagonal threshold", 1e-10);
489 printf(
"CMS: Reparing diagonal after next level generation\n");
493 OldEdgeLevel->
Set(
"NodeMatrix", OldSmootherMatrix);
495 OldSmootherMatrix = NewKn;
498 EdgeL->
Set(
"NodeMatrix", NodeAggMatrix);
502 EdgeL->
Set(
"NodeMatrix", NodeOp);
503 EdgeL->
Set(
"NodeAggMatrix", NodeOp);
506 OldEdgeLevel = EdgeL;
513 EdgeL->
Set(
"NodeImporter", importer);
519 std::string fine_smoother = precList11_.get<std::string>(
"smoother: type");
521 SM_Matrix_->setObjectLabel(
"A(1,1)");
522 precList11_.set(
"coarse: max size", 1);
523 precList11_.set(
"max levels", Hierarchy22_->GetNumLevels());
525 const bool refmaxwellCoarseSolve = (precList11_.get<std::string>(
"coarse: type",
526 MasterList::getDefault<std::string>(
"coarse: type")) ==
"RefMaxwell");
527 if (refmaxwellCoarseSolve) {
528 GetOStream(
Runtime0) <<
"Maxwell1::compute(): Will set up RefMaxwell coarse solver" << std::endl;
529 precList11_.set(
"coarse: type",
"none");
536 Hierarchy11_->setlib(SM_Matrix_->getDomainMap()->lib());
537 Hierarchy11_->SetProcRankVerbose(SM_Matrix_->getDomainMap()->getComm()->getRank());
540 mueLuFactory->SetupHierarchy(*Hierarchy11_);
542 if (refmaxwellCoarseSolve) {
543 GetOStream(
Runtime0) <<
"Maxwell1::compute(): Setting up RefMaxwell coarse solver" << std::endl;
544 auto coarseLvl = Hierarchy11_->GetLevel(Hierarchy11_->GetNumLevels() - 1);
546 precList11_.sublist(
"coarse: params")));
547 coarseSolver->Setup(*coarseLvl);
548 coarseLvl->Set(
"PreSmoother",
549 rcp_dynamic_cast<SmootherBase>(coarseSolver,
true));
552 if (mode_ == MODE_REFMAXWELL) {
553 if (Hierarchy11_->GetNumLevels() > 1) {
566 #ifdef MUELU_MAXWELL1_DEBUG
567 for (
int i = 0; i < Hierarchy11_->GetNumLevels(); i++) {
574 auto nrmE = EdgeMatrix->getFrobeniusNorm();
575 auto nrmN = NodeMatrix->getFrobeniusNorm();
576 auto nrmNa = NodeAggMatrix->getFrobeniusNorm();
577 auto nrmD0 = D0->getFrobeniusNorm();
579 std::cout <<
"DEBUG: Norms on Level " << i <<
" E/N/NA/D0 = " << nrmE <<
" / " << nrmN <<
" / " << nrmNa <<
" / " << nrmD0 << std::endl;
580 std::cout <<
"DEBUG: NNZ on Level " << i <<
" E/N/NA/D0 = " << EdgeMatrix->getGlobalNumEntries() <<
" / " << NodeMatrix->getGlobalNumEntries() <<
" / " << NodeAggMatrix->getGlobalNumEntries() <<
" / " << D0->getGlobalNumEntries() << std::endl;
585 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
589 RAPlist.
set(
"rap: fix zero diagonals",
false);
591 std::string labelstr =
"NodeMatrix (Level 0)";
596 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
598 if (mode_ == MODE_REFMAXWELL) {
601 residualFine_ = MultiVectorFactory::Build(SM_Matrix_->getRangeMap(), numVectors);
603 if (!Hierarchy11_.is_null() && Hierarchy11_->GetNumLevels() > 1) {
606 residual11c_ = MultiVectorFactory::Build(A->getRangeMap(), numVectors);
607 update11c_ = MultiVectorFactory::Build(A->getDomainMap(), numVectors);
610 if (!Hierarchy22_.is_null()) {
611 residual22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
612 update22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
617 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
619 if (dump_matrices_) {
620 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
625 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
627 if (dump_matrices_) {
628 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
633 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
635 if (dump_matrices_) {
636 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
641 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
643 if (dump_matrices_) {
644 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
645 std::ofstream out(name);
646 for (
size_t i = 0; i < Teuchos::as<size_t>(v.
size()); i++)
651 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
653 if (dump_matrices_) {
654 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
655 std::ofstream out(name);
656 auto vH = Kokkos::create_mirror_view(v);
657 Kokkos::deep_copy(vH, v);
658 for (
size_t i = 0; i < vH.size(); i++)
659 out << vH[i] <<
"\n";
663 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
675 return Teuchos::null;
678 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
680 bool reuse = !SM_Matrix_.is_null();
681 SM_Matrix_ = SM_Matrix_new;
682 dump(*SM_Matrix_,
"SM.m");
687 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
691 if (!allEdgesBoundary_ && X.
getNumVectors() != residualFine_->getNumVectors())
698 if (Fine->IsAvailable(
"PreSmoother")) {
701 preSmoo->
Apply(X, RHS,
true);
711 if (!P11_.is_null()) {
714 Hierarchy11_->Iterate(*residual11c_, *update11c_,
true, 1);
718 if (!allNodesBoundary_) {
721 Hierarchy22_->Iterate(*residual22_, *update22_,
true, 0);
729 if (!allNodesBoundary_)
734 if (Fine->IsAvailable(
"PostSmoother")) {
737 postSmoo->
Apply(X, RHS,
false);
741 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
743 Hierarchy11_->Iterate(RHS, X, 1,
true);
746 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
752 if (mode_ == MODE_GMHD_STANDARD)
753 HierarchyGmhd_->Iterate(RHS, X, 1,
true);
754 else if (mode_ == MODE_STANDARD || mode_ == MODE_EDGE_ONLY)
755 applyInverseStandard(RHS, X);
756 else if (mode_ == MODE_REFMAXWELL)
757 applyInverseRefMaxwellAdditive(RHS, X);
762 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
767 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
777 #ifdef HAVE_MUELU_DEBUG
779 TEUCHOS_ASSERT(Kn_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getDomainMap()));
780 TEUCHOS_ASSERT(Kn_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getDomainMap()));
783 TEUCHOS_ASSERT(D0_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getRowMap()));
786 Hierarchy11_ = Teuchos::null;
787 Hierarchy22_ = Teuchos::null;
788 HierarchyGmhd_ = Teuchos::null;
789 if (mode_ != MODE_GMHD_STANDARD) mode_ = MODE_STANDARD;
793 allEdgesBoundary_ =
false;
794 allNodesBoundary_ =
false;
795 dump_matrices_ =
false;
796 enable_reuse_ =
false;
798 applyBCsTo22_ =
false;
808 RCP<Matrix> D0copy = MatrixFactory::Build(D0_Matrix->getRowMap(), D0_Matrix->getColMap(), 0);
809 RCP<CrsMatrix> D0copyCrs = rcp_dynamic_cast<CrsMatrixWrap>(D0copy,
true)->getCrsMatrix();
813 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getAllValues(D0rowptr_RCP, D0colind_RCP, D0vals_RCP);
818 D0copyCrs->allocateAllValues(D0vals_RCP.
size(), D0copyrowptr_RCP, D0copycolind_RCP, D0copyvals_RCP);
819 D0copyrowptr_RCP.
deepCopy(D0rowptr_RCP());
820 D0copycolind_RCP.deepCopy(D0colind_RCP());
821 D0copyvals_RCP.
deepCopy(D0vals_RCP());
822 D0copyCrs->setAllValues(D0copyrowptr_RCP,
825 D0copyCrs->expertStaticFillComplete(D0_Matrix->getDomainMap(), D0_Matrix->getRangeMap(),
826 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getImporter(),
827 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getExporter());
830 D0_Matrix_ = MatrixFactory::BuildCopy(D0_Matrix);
832 Kn_Matrix_ = Kn_Matrix;
834 Nullspace_ = Nullspace;
836 if (!Kn_Matrix_.is_null()) dump(*Kn_Matrix_,
"Kn.m");
837 if (!Nullspace_.is_null()) dump(*Nullspace_,
"nullspace.m");
838 if (!Coords_.is_null()) dumpCoords(*Coords_,
"coords.m");
841 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
844 std::ostringstream oss;
850 if (!Kn_Matrix_.is_null())
851 root = comm->getRank();
860 oss <<
"\n--------------------------------------------------------------------------------\n"
861 <<
"--- Maxwell1 Summary ---\n"
862 "--------------------------------------------------------------------------------"
869 SM_Matrix_->getRowMap()->getComm()->barrier();
871 numRows = SM_Matrix_->getGlobalNumRows();
872 nnz = SM_Matrix_->getGlobalNumEntries();
887 oss <<
"block " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " << std::setw(9) <<
" nnz/row" << std::endl;
888 oss <<
"(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
890 if (!Kn_Matrix_.is_null()) {
892 numRows = Kn_Matrix_->getGlobalNumRows();
893 nnz = Kn_Matrix_->getGlobalNumEntries();
895 oss <<
"(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
900 std::string outstr = oss.str();
904 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
906 int strLength = outstr.size();
907 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
908 if (comm->getRank() != root)
909 outstr.resize(strLength);
910 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
915 if (!Hierarchy11_.is_null())
916 Hierarchy11_->describe(out, GetVerbLevel());
918 if (!Hierarchy22_.is_null())
919 Hierarchy22_->describe(out, GetVerbLevel());
921 if (!HierarchyGmhd_.is_null())
922 HierarchyGmhd_->describe(out, GetVerbLevel());
926 std::ostringstream oss2;
928 oss2 <<
"Sub-solver distribution over ranks" << std::endl;
929 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;
931 int numProcs = comm->getSize();
939 if (!Kn_Matrix_.is_null())
941 std::vector<char> states(numProcs, 0);
943 MPI_Gather(&status, 1, MPI_CHAR, &states[0], 1, MPI_CHAR, 0, *rawMpiComm);
945 states.push_back(status);
948 int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
949 for (
int proc = 0; proc < numProcs; proc += rowWidth) {
950 for (
int j = 0; j < rowWidth; j++)
951 if (proc + j < numProcs)
952 if (states[proc + j] == 0)
954 else if (states[proc + j] == 1)
956 else if (states[proc + j] == 2)
963 oss2 <<
" " << proc <<
":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
972 #define MUELU_MAXWELL1_SHORT
973 #endif // ifdef MUELU_MAXWELL1_DEF_HPP
Important warning messages (one line)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > PtAPWrapper(const RCP< Matrix > &A, const RCP< Matrix > &P, Teuchos::ParameterList ¶ms, std::string &label)
Performs an P^T AP.
static void ZeroDirichletCols(Teuchos::RCP< Matrix > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
const Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
virtual std::string getObjectLabel() const
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
void dump(const Matrix &A, std::string name) const
dump out matrix
static magnitudeType eps()
static void AddNonSerializableDataToHierarchy(HierarchyManager &HM, Hierarchy &H, const ParameterList &nonSerialList)
Add non-serializable data to Hierarchy.
T & get(const std::string &name, T def_value)
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
void GMHDSetupHierarchy(Teuchos::ParameterList &List) const
Sets up hiearchy for GMHD matrices that include generalized Ohms law equations.
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)
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
One-liner description of what is happening.
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::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReplaceNonZerosWithOnes(const RCP< Matrix > &original)
Creates a copy of a matrix where the non-zero entries are replaced by ones.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
static void CheckRepairMainDiagonal(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &Ac, bool const &repairZeroDiagonals, Teuchos::FancyOStream &fos, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType threshold=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< Scalar >::magnitudeType >::zero())
MsgType toVerbLevel(const std::string &verbLevelStr)
Print even more statistics.
static RCP< Time > getNewTimer(const std::string &name)
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...
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
bool isParameter(const std::string &name) const
bool remove(std::string const &name, bool throwIfNotExists=true)
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
void setParameters(Teuchos::ParameterList &list)
Set parameters.
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
bool isSublist(const std::string &name) const
Teuchos::RCP< Teuchos::TimeMonitor > getTimer(std::string name, RCP< const Teuchos::Comm< int > > comm=Teuchos::null) const
get a (synced) timer
void deepCopy(const ArrayView< const T > &av)
Class that encapsulates Operator smoothers.
void applyInverseRefMaxwellAdditive(const MultiVector &RHS, MultiVector &X) const
apply RefMaxwell additive 2x2 style cycle
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::memory_space > &BCrowsKokkos, Kokkos::View< bool *, typename Node::device_type::memory_space > &BCcolsKokkos, Kokkos::View< bool *, typename Node::device_type::memory_space > &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.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Print all timing information.
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 applyInverseStandard(const MultiVector &RHS, MultiVector &X) const
apply standard Maxwell1 cycle
const 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.
bool isType(const std::string &name) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
virtual size_t getNumVectors() const =0
Teuchos::RCP< Matrix > generate_kn() const
Generates the Kn matrix.
static std::string translate(Teuchos::ParameterList ¶mList, const std::string &defaultVals="")
: Translate ML parameters to MueLu parameter XML string
void dumpCoords(const RealValuedMultiVector &X, std::string name) const
dump out real-valued multivector
virtual void Apply(MultiVector &x, const MultiVector &rhs, bool InitialGuessIsZero=false) const =0
Apply smoother.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
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...
static void CopyBetweenHierarchies(Hierarchy &fromHierarchy, Hierarchy &toHierarchy, const std::string fromLabel, const std::string toLabel, const std::string dataType)