From 05c5c81d9e13a66593471030e8fb5f206a92f151 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Thu, 18 Jul 2024 07:55:52 -0600 Subject: [PATCH] MueLu: Fixing Maxwell1 so it robustly works w/ repartitioning (#13079) * MueLu Maxwell1 fixes (with lots of debug) * MueLu: This actually works. Stick a fork in the sand * MueLu: Reerting to remove debug statements * MueLu: Reverting to remove debug statements * MueLu: Reverting to remove debug statements * MueLu: Reverting to remove debug statements * MueLu: Reverting to remove debug statements * MueLu: The small test problem seems to be working. Stick a fork in it * MueLu: Removing debug output * MueLu: Updating maxwell unit tests * MueLu: Adding accessor to Hierarchy manager * MueLu: Pushing HM cleanup to Maxwell1 * MueLu: Removing comments * MueLu: Removing comments and output * MueLu: Removing comments and output * MueLu: Sticking a for in something working * MueLu: Removing extra gunk * MueLu: Removing debugging info * MueLu: Initial alloc cleanup * MueLu: caglusa's suggestion * MueLu: CLANGgit push origin develop:csiefer-maxwell1-fix! * MueLu: Fixing off-by-one error --- packages/muelu/doc/UsersGuide/masterList.xml | 24 +++---- .../muelu/doc/UsersGuide/paramlist_hidden.tex | 8 +-- .../src/Interface/MueLu_HierarchyManager.hpp | 20 ++++-- .../MueLu_ParameterListInterpreter_def.hpp | 11 +-- .../src/MueCentral/MueLu_Hierarchy_decl.hpp | 3 +- .../muelu/src/MueCentral/MueLu_MasterList.cpp | 6 +- .../src/Operators/MueLu_Maxwell1_decl.hpp | 1 + .../src/Operators/MueLu_Maxwell1_def.hpp | 71 +++++++++++++------ .../Operators/MueLu_Maxwell_Utils_decl.hpp | 2 +- .../src/Operators/MueLu_Maxwell_Utils_def.hpp | 2 +- .../MueLu_RepartitionFactory_def.hpp | 8 +++ .../MueLu_ReitzingerPFactory_def.hpp | 59 ++++++++++++--- packages/muelu/test/maxwell/CMakeLists.txt | 1 + packages/muelu/test/maxwell/Maxwell3D.cpp | 15 +++- .../maxwell/Maxwell_Reitzinger_repart.xml | 8 ++- 15 files changed, 174 insertions(+), 65 deletions(-) diff --git a/packages/muelu/doc/UsersGuide/masterList.xml b/packages/muelu/doc/UsersGuide/masterList.xml index 123b365d7900..0a8a6b4f247f 100644 --- a/packages/muelu/doc/UsersGuide/masterList.xml +++ b/packages/muelu/doc/UsersGuide/masterList.xml @@ -1053,22 +1053,12 @@ string "{}" Keep a subset of the hierarchy data after the - hierarchy is constructed. Use comma-separated array in braces + hierarchy is constructed. This could be for disk I/O, or for + use in other preconditioners. Use comma-separated array in braces false not supported by ML - - save data - string - "{}" - Data to dump to disk after the - hierarchy is constructed. Use comma-separated array in braces - false - not supported by ML - - - print initial parameters bool @@ -1786,6 +1776,16 @@ Only used when tentative: calculate qr is set to false. false + + repartition: save importer + bool + false + Mark Importer as 'Final' to keep it on the hierarchy after setup is complete. + not supported by ML + false + + + diff --git a/packages/muelu/doc/UsersGuide/paramlist_hidden.tex b/packages/muelu/doc/UsersGuide/paramlist_hidden.tex index 0bb3e508c615..3cce34b667ca 100644 --- a/packages/muelu/doc/UsersGuide/paramlist_hidden.tex +++ b/packages/muelu/doc/UsersGuide/paramlist_hidden.tex @@ -236,10 +236,8 @@ \textit{colmap\_X\_level.mm}.} \cbb{keep data}{string}{"{}"}{Keep a subset of the hierarchy data after the - hierarchy is constructed. Use comma-separated array in braces} - -\cbb{save data}{string}{"{}"}{Data to dump to disk after the - hierarchy is constructed. Use comma-separated array in braces} + hierarchy is constructed. This could be for disk I/O, or for + use in other preconditioners. Use comma-separated array in braces} \cbb{print initial parameters}{bool}{true}{Print parameters provided for a hierarchy construction.} @@ -401,6 +399,8 @@ \cbb{repartition: use subcommunicators}{bool}{true}{Use subcommunicators on coarser levels.} +\cbb{repartition: save importer}{bool}{false}{Mark Importer as 'Final' to keep it on the hierarchy after setup is complete.} + \cbb{rap: relative diagonal floor}{Array(double)}{{}}{Will boost diagonals on A matrices to the given floor.} \cbb{rap: fix zero diagonals}{bool}{false}{Set zero diagonals on coarse grids to one.} diff --git a/packages/muelu/src/Interface/MueLu_HierarchyManager.hpp b/packages/muelu/src/Interface/MueLu_HierarchyManager.hpp index 031fa5d73ef2..cf23ef5769be 100644 --- a/packages/muelu/src/Interface/MueLu_HierarchyManager.hpp +++ b/packages/muelu/src/Interface/MueLu_HierarchyManager.hpp @@ -171,6 +171,7 @@ class HierarchyManager : public HierarchyFactory>::const_iterator it = keep_.find(i); if (it != keep_.end()) { @@ -221,9 +223,9 @@ class HierarchyManager : public HierarchyFactory> FactoryMap; @@ -311,7 +320,6 @@ class HierarchyManager : public HierarchyFactory aggregatesToPrint_; Teuchos::Array elementToNodeMapsToPrint_; - // Data we'll need to save, not necessarily print - Teuchos::Array dataToSave_; + // Data we'll need to keep, either to dump to disk or to use post-setup + Teuchos::Array dataToKeep_; // Matrices we'll need to print std::map> matricesToPrint_; diff --git a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp index 8b3c8ea64115..609514f0572a 100644 --- a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp +++ b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp @@ -256,11 +256,11 @@ void ParameterListInterpreter:: (void)MUELU_TEST_AND_SET_VAR(paramList, "debug: graph level", int, this->graphOutputLevel_); - // Generic data saving (this saves the data on all levels) - if (paramList.isParameter("save data")) - this->dataToSave_ = Teuchos::getArrayFromStringParameter(paramList, "save data"); + // Generic data keeping (this keeps the data on all levels) + if (paramList.isParameter("keep data")) + this->dataToKeep_ = Teuchos::getArrayFromStringParameter(paramList, "keep data"); - // Save level data + // Export level data if (paramList.isSublist("export data")) { ParameterList printList = paramList.sublist("export data"); @@ -1630,6 +1630,7 @@ void ParameterListInterpreter:: MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: print partition distribution", bool, repartParams); MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: remap parts", bool, repartParams); MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: remap num values", int, repartParams); + MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "repartition: save importer", bool, repartParams); repartFactory->SetParameterList(repartParams); repartFactory->SetFactory("A", manager.GetFactory("A")); repartFactory->SetFactory("number of partitions", manager.GetFactory("number of partitions")); @@ -2272,7 +2273,7 @@ void ParameterListInterpreter:: } if (hieraList.isParameter("number of vectors")) { - this->numDesiredLevel_ = hieraList.get("number of vectors"); + this->sizeOfMultiVectors_ = hieraList.get("number of vectors"); hieraList.remove("number of vectors"); } diff --git a/packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp b/packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp index 35c4016df946..1ce9a84862b2 100644 --- a/packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp +++ b/packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp @@ -135,8 +135,9 @@ class Hierarchy : public BaseClass { template friend class Hierarchy; - private: int LastLevelID() const { return Levels_.size() - 1; } + + private: void DumpCurrentGraph(int level) const; public: diff --git a/packages/muelu/src/MueCentral/MueLu_MasterList.cpp b/packages/muelu/src/MueCentral/MueLu_MasterList.cpp index e0e9f0139168..a7c8176126ca 100644 --- a/packages/muelu/src/MueCentral/MueLu_MasterList.cpp +++ b/packages/muelu/src/MueCentral/MueLu_MasterList.cpp @@ -238,7 +238,6 @@ namespace MueLu { "" "" "" - "" "" "" "" @@ -317,6 +316,7 @@ namespace MueLu { "" "" "" + "" "" "" "" @@ -741,8 +741,6 @@ namespace MueLu { ("keep data","keep data") - ("save data","save data") - ("ML print initial list","print initial parameters") ("print unused","print unused parameters") @@ -899,6 +897,8 @@ namespace MueLu { ("repartition: use subcommunicators","repartition: use subcommunicators") + ("repartition: save importer","repartition: save importer") + ("rap: relative diagonal floor","rap: relative diagonal floor") ("rap: fix zero diagonals","rap: fix zero diagonals") diff --git a/packages/muelu/src/Operators/MueLu_Maxwell1_decl.hpp b/packages/muelu/src/Operators/MueLu_Maxwell1_decl.hpp index 4eb99c4343e1..61f9db95c339 100644 --- a/packages/muelu/src/Operators/MueLu_Maxwell1_decl.hpp +++ b/packages/muelu/src/Operators/MueLu_Maxwell1_decl.hpp @@ -19,6 +19,7 @@ #include "MueLu_Level_fwd.hpp" #include "MueLu_Hierarchy_fwd.hpp" #include "MueLu_RAPFactory_fwd.hpp" +#include "MueLu_RebalanceAcFactory_fwd.hpp" #include "MueLu_PerfUtils_fwd.hpp" #include "MueLu_SmootherBase_fwd.hpp" diff --git a/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp b/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp index 56b0c9bedd12..e962c9f3682c 100644 --- a/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp +++ b/packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp @@ -27,6 +27,7 @@ #include "MueLu_Level.hpp" #include "MueLu_Hierarchy.hpp" #include "MueLu_RAPFactory.hpp" +#include "MueLu_RebalanceAcFactory.hpp" #include "MueLu_Monitor.hpp" #include "MueLu_PerfUtils.hpp" #include "MueLu_ParameterListInterpreter.hpp" @@ -111,6 +112,7 @@ void Maxwell1::setParameters(Teuchos: list = newList; } + std::string mode_string = list.get("maxwell1: mode", MasterList::getDefault("maxwell1: mode")); applyBCsTo22_ = list.get("maxwell1: apply BCs to 22", false); dump_matrices_ = list.get("maxwell1: dump matrices", MasterList::getDefault("maxwell1: dump matrices")); @@ -273,16 +275,20 @@ void Maxwell1::compute(bool reuse) { if (!Coords_.is_null()) precList22_.sublist("user data").set("Coordinates", Coords_); - /* Repartitioning *must* be in sync between hierarchies, but the - only thing we need to watch here is the subcomms, since ReitzingerPFactory - won't look at all the other stuff */ + /* Repartitioning *must* be in sync between hierarchies, which means + that we need to keep the importers in sync too */ if (precList22_.isParameter("repartition: enable")) { bool repartition = precList22_.get("repartition: enable"); precList11_.set("repartition: enable", repartition); - // If we're repartitioning (2,2), we need to rebalance for (1,1) to do the right thing - if (repartition) + // If we're repartitioning (2,2), we need to rebalance for (1,1) to do the right thing, + // as well as keep the importers + if (repartition) { + // FIXME: We should probably update rather than clobber + precList22_.set("repartition: save importer", true); + // precList22_.set("repartition: rebalance P and R", false); precList22_.set("repartition: rebalance P and R", true); + } if (precList22_.isParameter("repartition: use subcommunicators")) { precList11_.set("repartition: use subcommunicators", precList22_.get("repartition: use subcommunicators")); @@ -407,8 +413,8 @@ void Maxwell1::compute(bool reuse) { Hierarchy11_->AddNewLevel(); RCP NodeL = Hierarchy22_->GetLevel(i); RCP EdgeL = Hierarchy11_->GetLevel(i); - RCP NodeOp = NodeL->Get >("A"); - RCP NodeAggMatrix = rcp_dynamic_cast(NodeOp); + RCP NodeAggOp = NodeL->Get >("A"); + RCP NodeAggMatrix = rcp_dynamic_cast(NodeAggOp); std::string labelstr = FormattingHelper::getColonLabel(EdgeL->getObjectLabel()); if (i == 0) { @@ -429,9 +435,16 @@ void Maxwell1::compute(bool reuse) { TEUCHOS_TEST_FOR_EXCEPTION(NodalP_ones.is_null(), Exceptions::RuntimeError, "Applying ones to prolongator failed"); EdgeL->Set("Pnodal", NodalP_ones); + // Get the importer if we have one (for repartitioning) + RCP importer; + if (NodeL->IsAvailable("Importer")) { + importer = NodeL->Get >("Importer"); + EdgeL->Set("NodeImporter", importer); + } + // If we repartition a processor away, a RCP is stuck // on the level instead of an RCP - if (!NodeAggMatrix.is_null()) { + if (!OldSmootherMatrix.is_null() && !NodalP_ones.is_null()) { EdgeL->Set("NodeAggMatrix", NodeAggMatrix); if (!have_generated_Kn) { // The user gave us a Kn on the fine level, so we're using a seperate aggregation @@ -446,12 +459,33 @@ void Maxwell1::compute(bool reuse) { Teuchos::ParameterList RAPlist; RAPlist.set("rap: fix zero diagonals", false); RCP NewKn = Maxwell_Utils::PtAPWrapper(OldSmootherMatrix, NodalP_ones, RAPlist, labelstr); - EdgeL->Set("NodeMatrix", NewKn); + + // If there's an importer, we need to Rebalance the NewKn + if (!importer.is_null()) { + ParameterList rebAcParams; + rebAcParams.set("repartition: use subcommunicators", true); + rebAcParams.set("repartition: use subcommunicators in place", true); + auto newAfact = rcp(new RebalanceAcFactory()); + newAfact->SetParameterList(rebAcParams); + RCP InPlaceMap = NodeAggMatrix.is_null() ? Teuchos::null : NodeAggMatrix->getRowMap(); + + Level fLevel, cLevel; + cLevel.SetPreviousLevel(Teuchos::rcpFromRef(fLevel)); + cLevel.Set("InPlaceMap", InPlaceMap); + cLevel.Set("A", NewKn); + cLevel.Request("A", newAfact.get()); + newAfact->Build(fLevel, cLevel); + + NewKn = cLevel.Get >("A", newAfact.get()); + EdgeL->Set("NodeMatrix", NewKn); + } else { + EdgeL->Set("NodeMatrix", NewKn); + } // Fix the old one double thresh = parameterList_.get("maxwell1: nodal smoother fix zero diagonal threshold", 1e-10); if (thresh > 0.0) { - printf("CMS: Reparing diagonal after next level generation\n"); + GetOStream(Runtime0) << "Maxwell1::compute(): Fixing zero diagonal for nodal smoother matrix" << std::endl; Scalar replacement = Teuchos::ScalarTraits::one(); Xpetra::MatrixUtils::CheckRepairMainDiagonal(OldSmootherMatrix, true, GetOStream(Warnings1), thresh, replacement); } @@ -464,19 +498,13 @@ void Maxwell1::compute(bool reuse) { } } else { // We've partitioned things away. - EdgeL->Set("NodeMatrix", NodeOp); - EdgeL->Set("NodeAggMatrix", NodeOp); + EdgeL->Set("NodeMatrix", NodeAggOp); + EdgeL->Set("NodeAggMatrix", NodeAggOp); } OldEdgeLevel = EdgeL; } - // Get the importer if we have one (for repartitioning) - // This will get used in ReitzingerPFactory - if (NodeL->IsAvailable("Importer")) { - auto importer = NodeL->Get >("Importer"); - EdgeL->Set("NodeImporter", importer); - } } // end Hierarchy22 loop //////////////////////////////////////////////////////////////////////////////// @@ -500,8 +528,11 @@ void Maxwell1::compute(bool reuse) { RCP > mueLuFactory = rcp(new ParameterListInterpreter(processedPrecList11, SM_Matrix_->getDomainMap()->getComm())); Hierarchy11_->setlib(SM_Matrix_->getDomainMap()->lib()); Hierarchy11_->SetProcRankVerbose(SM_Matrix_->getDomainMap()->getComm()->getRank()); - // Stick the non-serializible data on the hierarchy. - HierarchyUtils::AddNonSerializableDataToHierarchy(*mueLuFactory, *Hierarchy11_, nonSerialList11); + + // Stick the non-serializible data on the hierarchy and do setup + if (nonSerialList11.numParams() > 0) { + HierarchyUtils::AddNonSerializableDataToHierarchy(*mueLuFactory, *Hierarchy11_, nonSerialList11); + } mueLuFactory->SetupHierarchy(*Hierarchy11_); if (refmaxwellCoarseSolve) { diff --git a/packages/muelu/src/Operators/MueLu_Maxwell_Utils_decl.hpp b/packages/muelu/src/Operators/MueLu_Maxwell_Utils_decl.hpp index b907b99c36c8..a9b37ccc41cf 100644 --- a/packages/muelu/src/Operators/MueLu_Maxwell_Utils_decl.hpp +++ b/packages/muelu/src/Operators/MueLu_Maxwell_Utils_decl.hpp @@ -97,7 +97,7 @@ class Maxwell_Utils : public VerboseObject { //! Performs an P^T AP static RCP > - PtAPWrapper(const RCP& A, const RCP& P, Teuchos::ParameterList& params, std::string& label); + PtAPWrapper(const RCP& A, const RCP& P, Teuchos::ParameterList& params, const std::string& label); }; } // namespace MueLu diff --git a/packages/muelu/src/Operators/MueLu_Maxwell_Utils_def.hpp b/packages/muelu/src/Operators/MueLu_Maxwell_Utils_def.hpp index 58535a8aece2..5107553005a6 100644 --- a/packages/muelu/src/Operators/MueLu_Maxwell_Utils_def.hpp +++ b/packages/muelu/src/Operators/MueLu_Maxwell_Utils_def.hpp @@ -252,7 +252,7 @@ void Maxwell_Utils:: template RCP > Maxwell_Utils:: - PtAPWrapper(const RCP& A, const RCP& P, ParameterList& params, std::string& label) { + PtAPWrapper(const RCP& A, const RCP& P, ParameterList& params, const std::string& label) { Level fineLevel, coarseLevel; fineLevel.SetFactoryManager(null); coarseLevel.SetFactoryManager(null); diff --git a/packages/muelu/src/Rebalancing/MueLu_RepartitionFactory_def.hpp b/packages/muelu/src/Rebalancing/MueLu_RepartitionFactory_def.hpp index 855efd3f0135..d7a47734b7d9 100644 --- a/packages/muelu/src/Rebalancing/MueLu_RepartitionFactory_def.hpp +++ b/packages/muelu/src/Rebalancing/MueLu_RepartitionFactory_def.hpp @@ -53,6 +53,7 @@ RCP RepartitionFactoryset >("A", Teuchos::null, "Factory of the matrix A"); @@ -366,6 +367,13 @@ void RepartitionFactory::Build(Level& Set(currentLevel, "Importer", rowMapImporter); + // Importer saving + bool save_importer = pL.get("repartition: save importer"); + if (save_importer) { + currentLevel.Set("Importer", rowMapImporter, NoFactory::get()); + currentLevel.AddKeepFlag("Importer", NoFactory::get(), MueLu::Final); + currentLevel.RemoveKeepFlag("Importer", NoFactory::get(), MueLu::UserData); // FIXME: This is a hack + } // ====================================================================================================== // Print some data // ====================================================================================================== diff --git a/packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_ReitzingerPFactory_def.hpp b/packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_ReitzingerPFactory_def.hpp index b78beb8831b6..dbc354bb7f0e 100644 --- a/packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_ReitzingerPFactory_def.hpp +++ b/packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_ReitzingerPFactory_def.hpp @@ -179,8 +179,8 @@ void ReitzingerPFactory::BuildP(Level size_t Ne = EdgeMatrix->getLocalNumRows(); size_t Nn = NodeMatrix->getLocalNumRows(); - // Upper bound on local number of coarse edges - size_t max_edges = (NodeMatrix->getLocalNumEntries() + Nn + 1) / 2; + // Notinal upper bound on local number of coarse edges + size_t max_edges = PnT_D0T->getLocalNumEntries(); ArrayRCP D0_rowptr(Ne + 1); ArrayRCP D0_colind(max_edges); ArrayRCP D0_values(max_edges); @@ -261,6 +261,16 @@ void ReitzingerPFactory::BuildP(Level continue; } + // Exponential memory reallocation, if needed + if (current + 1 >= max_edges) { + max_edges *= 2; + D0_colind.resize(max_edges); + D0_values.resize(max_edges); + } + if (current / 2 + 1 >= D0_rowptr.size()) { + D0_rowptr.resize(2 * D0_rowptr.size() + 1); + } + // NOTE: "i" here might not be a valid local column id, so we read it from the map D0_colind[current] = local_column_i; D0_values[current] = -1; @@ -278,6 +288,11 @@ void ReitzingerPFactory::BuildP(Level D0_colind.resize(current); D0_values.resize(current); + // Handle empty ranks gracefully + if (num_coarse_edges == 0) { + D0_rowptr[0] = 0; + } + // Count the total number of edges // NOTE: Since we solve the ownership issue above, this should do what we want RCP ownedCoarseEdgeMap = Xpetra::MapFactory::Build(EdgeMatrix->getRowMap()->lib(), GO_INVALID, num_coarse_edges, EdgeMatrix->getRowMap()->getIndexBase(), EdgeMatrix->getRowMap()->getComm()); @@ -304,30 +319,31 @@ void ReitzingerPFactory::BuildP(Level #if 0 { char fname[80]; - printf("[%d] D0: ia.size() = %d ja.size() = %d\n",MyPID,(int)ia.size(),(int)ja.size()); - printf("[%d] D0: ia :",MyPID); + int fine_level_id = fineLevel.GetLevelID(); + printf("[%d] Level %d D0: ia.size() = %d ja.size() = %d\n",MyPID,fine_level_id,(int)ia.size(),(int)ja.size()); + printf("[%d] Level %d D0: ia :",MyPID,fine_level_id); for(int i=0; i<(int)ia.size(); i++) printf("%d ",(int)ia[i]); - printf("\n[%d] D0: global ja :",MyPID); + printf("\n[%d] Level %d D0: global ja :",MyPID,fine_level_id); for(int i=0; i<(int)ja.size(); i++) printf("%d ",(int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i])); - printf("\n[%d] D0: local ja :",MyPID); + printf("\n[%d] Level %d D0: local ja :",MyPID,fine_level_id); for(int i=0; i<(int)ja.size(); i++) printf("%d ",(int)ja[i]); printf("\n"); - sprintf(fname,"D0_global_ja_%d_%d.dat",MyPID,fineLevel.GetLevelID()); + sprintf(fname,"D0_global_ja_%d_%d.dat",MyPID,fine_level_id); FILE * f = fopen(fname,"w"); for(int i=0; i<(int)ja.size(); i++) fprintf(f,"%d ",(int)ownedPlusSharedCoarseNodeMap->getGlobalElement(ja[i])); fclose(f); - sprintf(fname,"D0_local_ja_%d_%d.dat",MyPID,fineLevel.GetLevelID()); + sprintf(fname,"D0_local_ja_%d_%d.dat",MyPID,fine_level_id); f = fopen(fname,"w"); for(int i=0; i<(int)ja.size(); i++) fprintf(f,"%d ",(int)ja[i]); fclose(f); - + } #endif D0_coarse->expertStaticFillComplete(ownedCoarseNodeMap, ownedCoarseEdgeMap); @@ -344,7 +360,28 @@ void ReitzingerPFactory::BuildP(Level RCP Pe; { SubFactoryMonitor m2(*this, "Generate Pe (pre-fix)", coarseLevel); - +#if 0 + { + // If you're concerned about processor / rank mismatches, this debugging code might help + int rank = D0->getRowMap()->getComm()->getRank(); + int fine_level = fineLevel.GetLevelID(); + printf("[%d] Level %d Checkpoint #2 Pn = %d/%d/%d/%d D0c = %d/%d/%d/%d D0 = %d/%d/%d/%d\n",rank,fine_level, + Pn->getRangeMap()->getComm()->getSize(), + Pn->getRowMap()->getComm()->getSize(), + Pn->getColMap()->getComm()->getSize(), + Pn->getDomainMap()->getComm()->getSize(), + D0_coarse_m->getRangeMap()->getComm()->getSize(), + D0_coarse_m->getRowMap()->getComm()->getSize(), + D0_coarse_m->getColMap()->getComm()->getSize(), + D0_coarse_m->getDomainMap()->getComm()->getSize(), + D0->getRangeMap()->getComm()->getSize(), + D0->getRowMap()->getComm()->getSize(), + D0->getColMap()->getComm()->getSize(), + D0->getDomainMap()->getComm()->getSize()); + fflush(stdout); + D0->getRowMap()->getComm()->barrier(); + } +#endif RCP dummy; RCP Pn_D0cT = XMM::Multiply(*Pn, false, *D0_coarse_m, true, dummy, out0, true, true, "Pn*D0c'", mm_params); @@ -422,6 +459,8 @@ void ReitzingerPFactory::BuildP(Level Set(coarseLevel, "Ptent", Pe); Set(coarseLevel, "D0", D0_coarse_m); + + /* This needs to be kept for the smoothers */ coarseLevel.Set("D0", D0_coarse_m, NoFactory::get()); coarseLevel.AddKeepFlag("D0", NoFactory::get(), MueLu::Final); coarseLevel.RemoveKeepFlag("D0", NoFactory::get(), MueLu::UserData); diff --git a/packages/muelu/test/maxwell/CMakeLists.txt b/packages/muelu/test/maxwell/CMakeLists.txt index 131ec81e3d91..9e70bcd24624 100644 --- a/packages/muelu/test/maxwell/CMakeLists.txt +++ b/packages/muelu/test/maxwell/CMakeLists.txt @@ -39,6 +39,7 @@ IF (${PACKAGE_NAME}_ENABLE_Belos) ARGS "--linAlgebra=Tpetra --reuse --xml=Maxwell.xml" ARGS "--linAlgebra=Tpetra --precType=MueLu-Maxwell1 --xml=Maxwell_Reitzinger.xml" ARGS "--linAlgebra=Tpetra --precType=MueLu-Maxwell1 --xml=Maxwell_Reitzinger_repart.xml" + ARGS "--linAlgebra=Tpetra --precType=MueLu-Maxwell1 --xml=Maxwell_Reitzinger_repart.xml --ensure-kn" COMM serial mpi NUM_MPI_PROCS 4 ) diff --git a/packages/muelu/test/maxwell/Maxwell3D.cpp b/packages/muelu/test/maxwell/Maxwell3D.cpp index 257e45454d42..e5929084cb25 100644 --- a/packages/muelu/test/maxwell/Maxwell3D.cpp +++ b/packages/muelu/test/maxwell/Maxwell3D.cpp @@ -29,6 +29,7 @@ // MueLu #include #include +#include #include #include @@ -512,6 +513,9 @@ int main_(Teuchos::CommandLineProcessor& clp, Xpetra::UnderlyingLib lib, int arg bool use_stacked_timer = false; clp.setOption("stacked-timer", "no-stacked-timer", &use_stacked_timer, "use stacked timer"); + bool ensure_kn = false; + clp.setOption("ensure-kn", + "no-ensure-kn", &ensure_kn, "generate a kn matrix if the user doesn't provide one"); std::string S_file, D0_file, M1_file, M0_file; if (!Teuchos::ScalarTraits::isComplex) { @@ -645,8 +649,17 @@ int main_(Teuchos::CommandLineProcessor& clp, Xpetra::UnderlyingLib lib, int arg Ms_Matrix = Xpetra::IO::Read(Ms_file, edge_map); else Ms_Matrix = M1_Matrix; - if (Kn_file != "") + if (Kn_file != "") { + *out << "User provided Kn matrix" << std::endl; Kn_Matrix = Xpetra::IO::Read(Kn_file, node_map); + } else if (ensure_kn) { + // The user requested we generate Kn_Matrix with an SpGEMM + Teuchos::ParameterList params; + *out << "Making Kn Matrix as requested by yser" << std::endl; + Kn_Matrix = MueLu::Maxwell_Utils::PtAPWrapper(SM_Matrix, D0_Matrix, params, "User Kn"); + } else { + *out << "NOT using a Kn matrix" << std::endl; + } if ((M0inv_file == "") && (M0_file != "")) { // nodal mass matrix diff --git a/packages/muelu/test/maxwell/Maxwell_Reitzinger_repart.xml b/packages/muelu/test/maxwell/Maxwell_Reitzinger_repart.xml index f5fd1cd49006..8266f50bf61f 100644 --- a/packages/muelu/test/maxwell/Maxwell_Reitzinger_repart.xml +++ b/packages/muelu/test/maxwell/Maxwell_Reitzinger_repart.xml @@ -40,6 +40,13 @@ + + + + + + + @@ -59,7 +66,6 @@ -