Merge pull request #7194 from mayrmt/7164-muelu-interface-node-mappin…

…g-user-data

MueLu: pass primal/dual interface node mapping via user data sublist

packages/muelu/src/Misc/MueLu_InterfaceAggregationFactory_def.hpp

@@ -61,22 +61,12 @@ RCP<const ParameterList> InterfaceAggregationFactory<Scalar, LocalOrdinal, Globa
 {
   RCP<ParameterList> validParamList = rcp(new ParameterList());
 
-  validParamList->set<RCP<const FactoryBase>>("A", null, "Generating factory of A");
-  validParamList->set<RCP<const FactoryBase>>("Aggregates", null, "Generating factory of the Aggregates (for block 0,0)");
-  validParamList->set<RCP<const FactoryBase>>("DualNodeID2PrimalNodeID", null, "Generating factory of the DualNodeID2PrimalNodeID map");
+  validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory of A");
+  validParamList->set<RCP<const FactoryBase>>("Aggregates", Teuchos::null, "Generating factory of the Aggregates (for block 0,0)");
+  validParamList->set<RCP<const FactoryBase>>("DualNodeID2PrimalNodeID", Teuchos::null, "Generating factory of the DualNodeID2PrimalNodeID map");
 
   validParamList->set<LocalOrdinal>("number of DOFs per dual node", Teuchos::ScalarTraits<LocalOrdinal>::one(), "Number of DOFs per dual node");
 
-  /*
-    DualNodeID2PrimalNodeID represents the mapping between the Dual Node IDs to the Primal Node IDs.
-    This map makes one assumption: the number of dofs per dual node is a constant and has to be set in
-    the parameter list as "number of DOFs per dual node".
-
-    Layout of the std::map<key, value>:
-    - key: local ID of dual node
-    - value: local ID of primal node
-  */
-  validParamList->set<RCP<std::map<LocalOrdinal, LocalOrdinal>>>("DualNodeID2PrimalNodeID - level 0", null, "");
   return validParamList;
 }
 
@@ -86,7 +76,17 @@ void InterfaceAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Dec
   Input(currentLevel, "A");
   Input(currentLevel, "Aggregates");
 
-  if (currentLevel.GetLevelID() != 0)
+  if (currentLevel.GetLevelID() == 0)
+  {
+    if(currentLevel.IsAvailable("DualNodeID2PrimalNodeID", NoFactory::get())) {
+      currentLevel.DeclareInput("DualNodeID2PrimalNodeID", NoFactory::get(), this);
+    } else {
+      TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("DualNodeID2PrimalNodeID", NoFactory::get()),
+                                 Exceptions::RuntimeError,
+                                 "DualNodeID2PrimalNodeID was not provided by the user on level 0!");
+    }
+  }
+  else
   {
     Input(currentLevel, "DualNodeID2PrimalNodeID");
   }
@@ -114,15 +114,9 @@ void InterfaceAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Bui
   RCP<Dual2Primal_type> lagr2dof;
 
   if (currentLevel.GetLevelID() == 0)
-  {
-    lagr2dof = pL.get<RCP<Dual2Primal_type>>("DualNodeID2PrimalNodeID - level 0");
-
-    TEUCHOS_TEST_FOR_EXCEPTION(lagr2dof.is_null(), std::runtime_error, prefix << " MueLu requires a provided DualNodeID2PrimalNodeID map on the finest level.");
-  }
+    lagr2dof = currentLevel.Get<RCP<Dual2Primal_type>>("DualNodeID2PrimalNodeID", NoFactory::get());
   else
-  {
     lagr2dof = Get<RCP<Dual2Primal_type>>(currentLevel, "DualNodeID2PrimalNodeID");
-  }
 
   const LocalOrdinal nDOFPerDualNode = pL.get<LocalOrdinal>("number of DOFs per dual node");
 

packages/muelu/src/MueCentral/MueLu_HierarchyUtils_decl.hpp

@@ -105,10 +105,30 @@ namespace MueLu {
 #include "MueLu_UseShortNames.hpp"
   public:
     /*!
-     * This routine is used by the CreateE/TpetraPreconditioner routines.
-     * Adds the following non-serializable data (A,P,R,Nullspace,Coordinates) from level-specific sublist nonSerialList,
-     * calling AddNewLevel as appropriate.
-     */
+    \brief Add non-serializable data to Hierarchy
+
+    Add non-serializable data given level-specific sublist \c nonSerialList to the Hierarchy \c H.
+    Calling \c AddLevel() along the way, if necessary.
+
+    Non-serializable data to be added:
+    - Operator "A"
+    - Prolongator "P"
+    - Restrictor "R"
+    - "M"
+    - "Mdiag"
+    - "K"
+    - Nullspace information "Nullspace"
+    - Coordinate information "Coordinates"
+    - "Node Comm"
+    - Primal-to-dual node mapping "DualNodeID2PrimalNodeID"
+    - "pcoarsen: element to node map
+
+    This routine is used by the CreateXpetraPreconditioner() routine.
+
+    @param HM
+    @param H
+    @param nonSerialList Parameter list containing non-serializable data
+    */
     static void AddNonSerializableDataToHierarchy(HierarchyManager& HM, Hierarchy& H, const ParameterList& nonSerialList);
   };
 

packages/muelu/src/MueCentral/MueLu_HierarchyUtils_def.hpp

@@ -94,7 +94,7 @@ namespace MueLu {
           const std::string& name = it2->first;
           TEUCHOS_TEST_FOR_EXCEPTION(name != "A" && name != "P" && name != "R" && name != "K"  && name != "M" && name != "Mdiag" &&
                                      name != "Nullspace" && name != "Coordinates" && name != "pcoarsen: element to node map" &&
-                                     name != "Node Comm" &&
+                                     name != "Node Comm" && name != "DualNodeID2PrimalNodeID" &&
                                      !IsParamMuemexVariable(name), Exceptions::InvalidArgument,
                                      std::string("MueLu::Utils::AddNonSerializableDataToHierarchy: parameter list contains unknown data type(") + name + ")");
 
@@ -131,7 +131,11 @@ namespace MueLu {
             level->AddKeepFlag(name,NoFactory::get(),MueLu::UserData);
             level->Set(name, Teuchos::getValue<RCP<const Teuchos::Comm<int> > >(it2->second), NoFactory::get());
           }
-
+          else if(name == "DualNodeID2PrimalNodeID")
+          {
+            level->AddKeepFlag(name,NoFactory::get(),MueLu::UserData);
+            level->Set(name, Teuchos::getValue<RCP<std::map<LO, LO>>>(it2->second), NoFactory::get());
+          }
 #ifdef HAVE_MUELU_INTREPID2
           else if (name == "pcoarsen: element to node map")
           {
@@ -187,7 +191,7 @@ namespace MueLu {
           const std::string& name = it2->first;
           TEUCHOS_TEST_FOR_EXCEPTION(name != "P" && name != "R"  && name != "K"  && name != "M" && name != "Mdiag" &&
                                      name != "Nullspace" && name != "Coordinates" && name != "pcoarsen: element to node map" && 
-                                     name != "Node Comm" &&
+                                     name != "Node Comm" && name != "DualNodeID2PrimalNodeID" &&
                                      !IsParamValidVariable(name), Exceptions::InvalidArgument,
                                      std::string("MueLu::Utils::AddNonSerializableDataToHierarchy: user data parameter list contains unknown data type (") + name + ")");
           if( name == "P" || name == "R" || name == "K" || name == "M") {
@@ -210,6 +214,11 @@ namespace MueLu {
             level->AddKeepFlag(name,NoFactory::get(),MueLu::UserData);
             level->Set(name, Teuchos::getValue<RCP<const Teuchos::Comm<int> > >(it2->second), NoFactory::get());
           }
+          else if(name == "DualNodeID2PrimalNodeID")
+          {
+            level->AddKeepFlag(name,NoFactory::get(),MueLu::UserData);
+            level->Set(name, Teuchos::getValue<RCP<std::map<LO, LO>>>(it2->second), NoFactory::get());
+          }
 #ifdef HAVE_MUELU_INTREPID2
           else if (name == "pcoarsen: element to node map")
           {

packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp

@@ -305,7 +305,7 @@ namespace MueLu {
     //@{
 
     //! Return a simple one-line description of this object.
-    std::string description() const; 
+    std::string description() const;
 
     /*! @brief Print the Hierarchy with some verbosity level to a FancyOStream object.
 
@@ -315,7 +315,7 @@ namespace MueLu {
     void describe(Teuchos::FancyOStream& out, const VerbLevel verbLevel = Default) const;
     void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_HIGH) const;
 
-    // Hierarchy::print is local hierarchy function, thus the statistics can be different from global ones
+    //! Hierarchy::print is local hierarchy function, thus the statistics can be different from global ones
     void print(std::ostream& out = std::cout, const VerbLevel verbLevel = (MueLu::Parameters | MueLu::Statistics0)) const;
 
     /*! Indicate whether the multigrid method is a preconditioner or a solver.
@@ -335,7 +335,7 @@ namespace MueLu {
     void setlib(Xpetra::UnderlyingLib inlib) { lib_ = inlib; }
     Xpetra::UnderlyingLib lib() { return lib_; }
 
-    // force recreation of cached description_ next time description() is called: 
+    //! force recreation of cached description_ next time description() is called:
     void ResetDescription() {
       description_ = "";
     }
@@ -355,59 +355,62 @@ namespace MueLu {
     //! Container for Level objects
     Array<RCP<Level> > Levels_;
 
-    // We replace coordinates GIDs to make them consistent with matrix GIDs,
-    // even if user does not do that.  Ideally, though, we should completely
-    // remove any notion of coordinate GIDs, and deal only with LIDs, assuming
-    // that they are consistent with matrix block IDs
+    //! We replace coordinates GIDs to make them consistent with matrix GIDs,
+    //! even if user does not do that.  Ideally, though, we should completely
+    //! remove any notion of coordinate GIDs, and deal only with LIDs, assuming
+    //! that they are consistent with matrix block IDs
     void ReplaceCoordinateMap(Level& level);
 
-    // Minimum size of a matrix on any level. If we fall below that, we stop
-    // the coarsening
+    //! Minimum size of a matrix on any level. If we fall below that, we stop
+    //! the coarsening
     Xpetra::global_size_t maxCoarseSize_;
 
-    // Potential speed up of the setup by skipping R construction, and using
-    // transpose matrix-matrix product for RAP
+    //! Potential speed up of the setup by skipping R construction, and using
+    //! transpose matrix-matrix product for RAP
     bool implicitTranspose_;
 
-    // Potential speed up of the solve by fusing prolongation and update steps.
-    // This can lead to more iterations to round-off error accumulation.
+    //! Potential speed up of the solve by fusing prolongation and update steps.
+    //! This can lead to more iterations to round-off error accumulation.
     bool fuseProlongationAndUpdate_;
 
-    // Potential speed up of the setup by skipping rebalancing of P and R, and
-    // doing extra import during solve
+    //! Potential speed up of the setup by skipping rebalancing of P and R, and
+    //! doing extra import during solve
     bool doPRrebalance_;
 
-    // Hierarchy may be used in a standalone mode, or as a preconditioner
+    //! Hierarchy may be used in a standalone mode, or as a preconditioner
     bool isPreconditioner_;
 
-    // V- or W-cycle
+    //! V- or W-cycle
     CycleType Cycle_;
 
-    // Scaling factor to be applied to coarse grid correction.
+    //! Scaling factor to be applied to coarse grid correction.
     double scalingFactor_;
 
-    // Epetra/Tpetra mode
+    //! Epetra/Tpetra mode
     Xpetra::UnderlyingLib lib_;
 
-    // cache description to avoid recreating in each call to description() - use ResetDescription() to force recreation in Setup, SetupRe, etc.
+    //! cache description to avoid recreating in each call to description() - use ResetDescription() to force recreation in Setup, SetupRe, etc.
     mutable std::string description_ = ""; // mutable so that we can lazily initialize in description(), which is declared const
 
-    //! Graph dumping
-    // If enabled, we dump the graph on a specified level into a specified file
+    /*!
+    @brief Graph dumping
+
+    If enabled, we dump the graph on a specified level into a specified file
+    */
     bool isDumpingEnabled_;
     int  dumpLevel_;
     std::string dumpFile_;
 
     //! Convergece rate
     MagnitudeType rate_;
 
-    // Level managers used during the Setup
+    //! Level managers used during the Setup
     Array<RCP<const FactoryManagerBase> > levelManagers_;
 
-    // Caching (Multi)Vectors used in Hierarchy::Iterate()
+    //! Caching (Multi)Vectors used in Hierarchy::Iterate()
     int sizeOfAllocatedLevelMultiVectors_;
     Array<RCP<MultiVector> > residual_, coarseRhs_, coarseX_, coarseImport_, coarseExport_, correction_;
-    
+
 
   }; //class Hierarchy
 

packages/muelu/src/Utils/MueLu_Utilities.cpp

@@ -62,12 +62,6 @@
 
 namespace MueLu {
 
-  /* Removes the following non-serializable data (A,P,R,Nullspace,Coordinates)
-     from level-specific sublists from inList
-     and moves it to nonSerialList.  Everything else is copied to serialList.
-     This function returns the level number of the highest level for which
-     non-serializable data was provided.
-  */
   long ExtractNonSerializableData(const Teuchos::ParameterList& inList, Teuchos::ParameterList& serialList, Teuchos::ParameterList& nonSerialList) {
     using Teuchos::ParameterList;
 
@@ -93,7 +87,7 @@ namespace MueLu {
         for (ParameterList::ConstIterator it2 = levelList.begin(); it2 != levelList.end(); it2++) {
           const std::string& name = it2->first;
           if (name == "A" || name == "P" || name == "R"  || name== "M" || name == "Mdiag" || name == "K" || name == "Nullspace" || name == "Coordinates"
-              || name == "Node Comm"
+              || name == "Node Comm" || name == "DualNodeID2PrimalNodeID"
 #ifdef HAVE_MUELU_INTREPID2 // For the IntrepidPCoarsenFactory
               || name == "pcoarsen: element to node map"
 #endif

packages/muelu/src/Utils/MueLu_Utilities_decl.hpp

@@ -895,10 +895,34 @@ namespace MueLu {
 
 
 
-  /*! Removes the following non-serializable data (A,P,R,Nullspace,Coordinates) from level-specific sublists from inList
-    and moves it to nonSerialList.  Everything else is copied to serialList.  This function returns the level number of the highest level
-    for which non-serializable data was provided.
-    */
+  /*!
+  \brief Extract non-serializable data from level-specific sublists and move it to a separate parameter list
+
+  Look through the level-specific sublists form \c inList, extract non-serializable data and move it to \c nonSerialList.
+  Everything else is copied to the \c serialList.
+
+  \note Data is considered "non-serializable" if it is not the same on every rank/processor.
+
+  Non-serializable data to be moved:
+  - Operator "A"
+  - Prolongator "P"
+  - Restrictor "R"
+  - "M"
+  - "Mdiag"
+  - "K"
+  - Nullspace information "Nullspace"
+  - Coordinate information "Coordinates"
+  - "Node Comm"
+  - Primal-to-dual node mapping "DualNodeID2PrimalNodeID"
+  - "pcoarsen: element to node map
+
+  @param[in] inList List with all input parameters/data as provided by the user
+  @param[out] serialList All serializable data from the input list
+  @param[out] nonSerialList All non-serializable, i.e. rank-specific data from the input list
+
+  @return This function returns the level number of the highest level for which non-serializable data was provided.
+
+  */
   long ExtractNonSerializableData(const Teuchos::ParameterList& inList, Teuchos::ParameterList& serialList, Teuchos::ParameterList& nonSerialList);
 
 

packages/muelu/test/meshtying/MeshTyingBlocked_SimpleSmoother.cpp

@@ -45,6 +45,7 @@
 // @HEADER
 
 // MueLu
+#include <MueLu_CreateXpetraPreconditioner.hpp>
 #include <MueLu_ConfigDefs.hpp>
 #include <MueLu_ParameterListInterpreter.hpp>
 
@@ -209,17 +210,11 @@ int main(int argc, char *argv[])
   std::string xmlFile = "myXML.xml";
 
   RCP<ParameterList> params = Teuchos::getParametersFromXmlFile(xmlFile);
-  RCP<ParameterList> paramsF, paramsI;
-  paramsF = sublist(params, "Factories");
-  paramsI = sublist(paramsF, "myInterfaceAggs2");
-  paramsI->set<RCP<std::map<LO, LO>>>("DualNodeID2PrimalNodeID - level 0", rcpFromRef(myLagr2Dof));
-
-  ParameterListInterpreter mueLuFactory(*params, comm);
-  RCP<Hierarchy> H = mueLuFactory.CreateHierarchy();
-  H->GetLevel(0)->Set("A", rcp_dynamic_cast<xpetra_matrix_type>(blockedMatrix));
+  ParameterList& userDataParams = params->sublist("user data");
+  userDataParams.set<RCP<std::map<LO, LO>>>("DualNodeID2PrimalNodeID", rcpFromRef(myLagr2Dof));
+
+  RCP<Hierarchy> H = MueLu::CreateXpetraPreconditioner(Teuchos::rcp_dynamic_cast<Matrix>(blockedMatrix, true), *params);
   H->IsPreconditioner(true);
-  H->SetDefaultVerbLevel(MueLu::Extreme);
-  mueLuFactory.SetupHierarchy(*H);
 
   // Create the preconditioned GMRES solver
   typedef xpetra_mvector_type MV;

packages/muelu/test/meshtying/myXML.xml

@@ -217,7 +217,7 @@
 
     <Parameter name="max levels"          type="int"      value="3"/>
     <Parameter name="coarse: max size"    type="int"      value="25"/>
-    <Parameter name="verbosity"           type="string"   value="High"/>
+    <Parameter name="verbosity"           type="string"   value="Extreme"/>
     <Parameter name="use kokkos refactor" type="bool"     value="false"/>
 
     <ParameterList name="AllLevel">