Loading [MathJax]/extensions/tex2jax.js
MueLu  Version of the Day
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_MultiPhys_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_MULTIPHYS_DEF_HPP
11 #define MUELU_MULTIPHYS_DEF_HPP
12 
13 #include <sstream>
14 #include "MueLu_ConfigDefs.hpp"
15 
16 #include "Xpetra_Map.hpp"
18 #include "Xpetra_MatrixUtils.hpp"
19 
20 #include "MueLu_MultiPhys_decl.hpp"
21 #include "MueLu_SaPFactory.hpp"
22 #include "MueLu_AggregationExportFactory.hpp"
23 #include "MueLu_Utilities.hpp"
24 #include "MueLu_Level.hpp"
25 #include "MueLu_Hierarchy.hpp"
26 #include "MueLu_RAPFactory.hpp"
27 #include "MueLu_Monitor.hpp"
28 #include "MueLu_PerfUtils.hpp"
29 #include "MueLu_ParameterListInterpreter.hpp"
31 #include <MueLu_HierarchyUtils.hpp>
32 #include "MueLu_VerbosityLevel.hpp"
35 
36 #ifdef HAVE_MUELU_CUDA
37 #include "cuda_profiler_api.h"
38 #endif
39 
40 namespace MueLu {
41 
42 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
44  return AmatMultiphysics_->getDomainMap();
45 }
46 
47 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
49  return AmatMultiphysics_->getRangeMap();
50 }
51 
52 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
54  // this operator only makes sense for the Combo when using TransP for R
55 
56  list.set("multigrid algorithm", "combine");
57  list.set("combine: numBlks", nBlks_);
58 
59  // Make sure verbosity gets passed to the sublists
60  std::string verbosity = list.get("verbosity", "high");
62 
63  arrayOfParamLists_.resize(nBlks_);
64  for (int i = 0; i < nBlks_; i++) {
65  std::string listName = "subblockList" + Teuchos::toString(i);
66  if (list.isSublist(listName)) {
67  arrayOfParamLists_[i] = Teuchos::rcpFromRef(list.sublist(listName));
68  } else
69  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Must provide sublist " + listName);
70 
71  arrayOfParamLists_[i]->set("verbosity", arrayOfParamLists_[i]->get("verbosity", verbosity));
72  if (OmitSubblockSmoother_) {
73  arrayOfParamLists_[i]->set("smoother: pre or post", "none");
74  arrayOfParamLists_[i]->set("smoother: type", "none");
75  }
76  }
77 
78  // Are we using Kokkos?
79  useKokkos_ = !Node::is_serial;
80  useKokkos_ = list.get("use kokkos refactor", useKokkos_);
81 
82  paramListMultiphysics_ = Teuchos::rcpFromRef(list);
83 }
84 
85 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
87  /*
88 
89  Create a set of AMG hierarchies whose interpolation matrices are used to build on combined
90  AMG hierarchy for a multiphysics problem
91 
92  */
93 
94  //#ifdef HAVE_MUELU_CUDA
95  // if (paramListMultiphysics_.get<bool>("multiphysics: cuda profile setup", false)) cudaProfilerStart();
96  //#endif
97 
98  std::string timerLabel;
99  if (reuse)
100  timerLabel = "MueLu MultiPhys: compute (reuse)";
101  else
102  timerLabel = "MueLu MultiPhys: compute";
103  RCP<Teuchos::TimeMonitor> tmCompute = getTimer(timerLabel);
104 
106  // Generate the (iii,iii) Hierarchy
107 
108  for (int iii = 0; iii < nBlks_; iii++) {
109  if (arrayOfCoords_ != Teuchos::null) {
110  if (arrayOfCoords_[iii] != Teuchos::null) {
111  arrayOfParamLists_[iii]->sublist("user data").set("Coordinates", arrayOfCoords_[iii]);
112  }
113  }
114 
115  if (arrayOfMaterials_ != Teuchos::null) {
116  if (arrayOfMaterials_[iii] != Teuchos::null) {
117  arrayOfParamLists_[iii]->sublist("user data").set("Material", arrayOfMaterials_[iii]);
118  }
119  }
120 
121  bool wantToRepartition = false;
122  if (paramListMultiphysics_->isParameter("repartition: enable"))
123  wantToRepartition = paramListMultiphysics_->get<bool>("repartition: enable");
124 
125  arrayOfParamLists_[iii]->set("repartition: enable", wantToRepartition);
126  arrayOfParamLists_[iii]->set("repartition: rebalance P and R", true);
127  arrayOfParamLists_[iii]->set("repartition: explicit via new copy rebalance P and R", true);
128 
129  if (paramListMultiphysics_->isParameter("repartition: use subcommunicators"))
130  arrayOfParamLists_[iii]->set("repartition: use subcommunicators", paramListMultiphysics_->isParameter("repartition: use subcommunicators"));
131  else
132  arrayOfParamLists_[iii]->set("repartition: use subcommunicators", true);
133  }
134  // repartitioning should only happen when createing the individual P's , not
135  // when combiing them
136 
137  paramListMultiphysics_->set<bool>("repartition: enable", false);
138 
139  LO maxLevels = 9999;
140  for (int i = 0; i < nBlks_; i++) {
141  std::string operatorLabel = "MultiPhys (" + Teuchos::toString(i) + "," + Teuchos::toString(i) + ")";
142  arrayOfAuxMatrices_[i]->setObjectLabel(operatorLabel);
143  arrayOfHierarchies_[i] = MueLu::CreateXpetraPreconditioner(arrayOfAuxMatrices_[i], *arrayOfParamLists_[i]);
144  LO tempNlevels = arrayOfHierarchies_[i]->GetGlobalNumLevels();
145  if (tempNlevels < maxLevels) maxLevels = tempNlevels;
146  }
147 
148  hierarchyMultiphysics_ = rcp(new Hierarchy("Combo"));
149  for (LO i = 0; i < maxLevels; i++) {
150  hierarchyMultiphysics_->AddNewLevel();
151  }
152  for (int i = 0; i < nBlks_; i++) {
153  std::string subblkName = "Psubblock" + Teuchos::toString(i);
154  MueLu::HierarchyUtils<SC, LO, GO, NO>::CopyBetweenHierarchies(*(arrayOfHierarchies_[i]), *(hierarchyMultiphysics_), "P", subblkName, "RCP<Matrix>");
155  }
156  paramListMultiphysics_->set("coarse: max size", 1);
157  paramListMultiphysics_->set("max levels", maxLevels);
158 
159  AmatMultiphysics_->setObjectLabel("A: block " + Teuchos::toString(nBlks_) + " x " + Teuchos::toString(nBlks_) + "multiphysics matrix");
160 
161  // Rip off non-serializable data before validation
162  Teuchos::ParameterList nonSerialListMultiphysics, processedListMultiphysics;
163  MueLu::ExtractNonSerializableData(*paramListMultiphysics_, processedListMultiphysics, nonSerialListMultiphysics);
164 
165  // Rip off the subblock List stuff as we don't need it any more and I think it messes up validator
166 
167  Teuchos::ParameterList stripped;
168  for (ParameterList::ConstIterator inListEntry = processedListMultiphysics.begin(); inListEntry != processedListMultiphysics.end(); inListEntry++) {
169  const std::string& levelName = inListEntry->first;
170  if (levelName.find("subblockList") != 0) stripped.setEntry(inListEntry->first, inListEntry->second);
171  }
172 
173  RCP<HierarchyManager<SC, LO, GO, NO>> mueLuFactory = rcp(new ParameterListInterpreter<SC, LO, GO, NO>(stripped, AmatMultiphysics_->getDomainMap()->getComm()));
174  hierarchyMultiphysics_->setlib(AmatMultiphysics_->getDomainMap()->lib());
175  hierarchyMultiphysics_->SetProcRankVerbose(AmatMultiphysics_->getDomainMap()->getComm()->getRank());
176 
177  // We don't need nullspace or coordinates, since we don't use them when just combining prolongators that have been already created
178  hierarchyMultiphysics_->GetLevel(0)->Set("A", AmatMultiphysics_);
179 
180  // Stick the non-serializible data on the hierarchy.
181  // Not sure that we need this, since we don't use it in building the multiphysics hierarchy
182  HierarchyUtils<SC, LO, GO, NO>::AddNonSerializableDataToHierarchy(*mueLuFactory, *hierarchyMultiphysics_, nonSerialListMultiphysics);
183  mueLuFactory->SetupHierarchy(*hierarchyMultiphysics_);
184 
185  describe(GetOStream(Runtime0));
186 }
187 
188 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
190  if (IsPrint(Timings)) {
191  if (!syncTimers_)
193  else {
194  if (comm.is_null())
195  return Teuchos::rcp(new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), AmatMultiphysics_->getRowMap()->getComm().ptr()));
196  else
198  }
199  } else
200  return Teuchos::null;
201 }
202 
203 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
205  bool reuse = !AmatMultiphysics_.is_null();
206  AmatMultiphysics_ = AmatMultiphysics_new;
207  if (ComputePrec) compute(reuse);
208 }
209 
210 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
212  hierarchyMultiphysics_->Iterate(RHS, X, 1, true);
213 }
214 
215 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
217  Teuchos::ETransp /* mode */,
218  Scalar /* alpha */,
219  Scalar /* beta */) const {
220  RCP<Teuchos::TimeMonitor> tm = getTimer("MueLu MultiPhys: solve");
221  hierarchyMultiphysics_->Iterate(RHS, X, 1, true);
222 }
223 
224 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
226  return false;
227 }
228 
229 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
231  initialize(const Teuchos::RCP<Matrix>& AmatMultiPhysics,
232  const Teuchos::ArrayRCP<RCP<Matrix>> arrayOfAuxMatrices,
233  const Teuchos::ArrayRCP<Teuchos::RCP<MultiVector>> arrayOfNullspaces,
235  const int nBlks,
237  const Teuchos::ArrayRCP<Teuchos::RCP<MultiVector>> arrayOfMaterials) {
238  arrayOfHierarchies_.resize(nBlks_);
239  for (int i = 0; i < nBlks_; i++) arrayOfHierarchies_[i] = Teuchos::null;
240 
241  // Default settings
242  useKokkos_ = false;
243  enable_reuse_ = false;
244  syncTimers_ = false;
245 
246  // set parameters
247  setParameters(List);
248 
249 } // initialize
250 
251 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
253  describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel /* verbLevel */) const {
254  std::ostringstream oss;
255 
256  RCP<const Teuchos::Comm<int>> comm = AmatMultiphysics_->getDomainMap()->getComm();
257 
258  oss << "\n--------------------------------------------------------------------------------\n"
259  << "--- MultiPhysics Summary ---\n"
260  "--------------------------------------------------------------------------------"
261  << std::endl;
262  oss << std::endl;
263 
264  GlobalOrdinal numRows;
265  GlobalOrdinal nnz;
266 
267  AmatMultiphysics_->getRowMap()->getComm()->barrier();
268 
269  for (int i = 0; i < nBlks_; i++) {
270  numRows = arrayOfAuxMatrices_[i]->getGlobalNumRows();
271  nnz = arrayOfAuxMatrices_[i]->getGlobalNumEntries();
272  Xpetra::global_size_t tt = numRows;
273  int rowspacer = 3;
274  while (tt != 0) {
275  tt /= 10;
276  rowspacer++;
277  }
278  tt = nnz;
279  int nnzspacer = 2;
280  while (tt != 0) {
281  tt /= 10;
282  nnzspacer++;
283  }
284 
285  oss << "block " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz " << std::setw(9) << " nnz/row" << std::endl;
286  oss << "(" << Teuchos::toString(i) << ", " << Teuchos::toString(i) << ")" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
287  }
288  oss << std::endl;
289 
290  out << oss.str();
291 
292  for (int i = 0; i < nBlks_; i++) {
293  arrayOfHierarchies_[i]->describe(out, GetVerbLevel());
294  }
295 
296 } // describe
297 
298 } // namespace MueLu
299 
300 #define MUELU_MULTIPHYS_SHORT
301 #endif // ifdef MUELU_MULTIPHYS_DEF_HPP
ParameterList & setEntry(const std::string &name, U &&entry)
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
ConstIterator end() const
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
LocalOrdinal LO
One-liner description of what is happening.
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
void setParameters(Teuchos::ParameterList &list)
Set parameters.
MsgType toVerbLevel(const std::string &verbLevelStr)
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...
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
void applyInverse(const MultiVector &RHS, MultiVector &X) const
apply standard MultiPhys cycle
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
const Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
bool isSublist(const std::string &name) const
params_t::ConstIterator ConstIterator
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
ConstIterator begin() const
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
size_t global_size_t
Print all timing information.
void initialize(const Teuchos::RCP< Matrix > &AmatMultiPhysics, const Teuchos::ArrayRCP< RCP< Matrix >> arrayOfAuxMatrices, const Teuchos::ArrayRCP< Teuchos::RCP< MultiVector >> arrayOfNullspaces, const Teuchos::ArrayRCP< Teuchos::RCP< RealValuedMultiVector >> arrayOfCoords, const int nBlks, Teuchos::ParameterList &List, const Teuchos::ArrayRCP< Teuchos::RCP< MultiVector >> arrayOfMaterials)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
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
Exception throws to report errors in the internal logical of the program.
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
void compute(bool reuse=false)
Setup the preconditioner.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
std::string toString(const T &t)
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)