MueLu  Version of the Day
MueLu_Maxwell1_def.hpp
Go to the documentation of this file.
1 //
2 // ***********************************************************************
3 //
4 // MueLu: A package for multigrid based preconditioning
5 // Copyright 2012 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact
38 // Jonathan Hu (jhu@sandia.gov)
39 // Andrey Prokopenko (aprokop@sandia.gov)
40 // Ray Tuminaro (rstumin@sandia.gov)
41 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 #ifndef MUELU_MAXWELL1_DEF_HPP
46 #define MUELU_MAXWELL1_DEF_HPP
47 
48 #include <sstream>
49 
50 #include "MueLu_ConfigDefs.hpp"
51 
52 #include "Xpetra_Map.hpp"
53 #include "Xpetra_MatrixMatrix.hpp"
54 #include "Xpetra_TripleMatrixMultiply.hpp"
55 #include "Xpetra_CrsMatrixUtils.hpp"
56 #include "Xpetra_MatrixUtils.hpp"
57 
58 #include "MueLu_Maxwell1_decl.hpp"
59 #include "MueLu_Maxwell_Utils.hpp"
60 
61 #include "MueLu_ReitzingerPFactory.hpp"
62 #include "MueLu_SaPFactory.hpp"
63 #include "MueLu_AggregationExportFactory.hpp"
64 #include "MueLu_Utilities.hpp"
65 #include "MueLu_Level.hpp"
66 #include "MueLu_Hierarchy.hpp"
67 #include "MueLu_RAPFactory.hpp"
68 #include "MueLu_PerfUtils.hpp"
69 #include "MueLu_ParameterListInterpreter.hpp"
71 
72 #if defined(HAVE_MUELU_KOKKOS_REFACTOR)
73 #include "MueLu_Utilities_kokkos.hpp"
74 #endif
75 
76 
77 #include "MueLu_VerbosityLevel.hpp"
78 
81 
82 #ifdef HAVE_MUELU_CUDA
83 #include "cuda_profiler_api.h"
84 #endif
85 
86 
87 namespace MueLu {
88 
89 
90  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
91  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getDomainMap() const {
92  return SM_Matrix_->getDomainMap();
93  }
94 
95 
96  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
97  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getRangeMap() const {
98  return SM_Matrix_->getRangeMap();
99  }
100 
101 
102  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
104 
105  std::string mode_string = list.get("maxwell1: mode", MasterList::getDefault<std::string>("maxwell1: mode"));
106  applyBCsTo22_ = list.get("maxwell1: apply BCs to 22", true);
107 
108  // Default smoother. We'll copy this around.
109  Teuchos::ParameterList defaultSmootherList;
110  defaultSmootherList.set("smoother: type", "CHEBYSHEV");
111  defaultSmootherList.sublist("smoother: params").set("chebyshev: degree",2);
112  defaultSmootherList.sublist("smoother: params").set("chebyshev: ratio eigenvalue",7.0);
113  defaultSmootherList.sublist("smoother: params").set("chebyshev: eigenvalue max iterations",30);
114 
115  // Make sure verbosity gets passed to the sublists
116  std::string verbosity = list.get("verbosity","high");
117 
118  // Check the validity of the run mode
119  if(mode_string == "standard") mode_ = MODE_STANDARD;
120  else if(mode_string == "refmaxwell") mode_ = MODE_REFMAXWELL;
121  else if(mode_string == "edge only") mode_ = MODE_EDGE_ONLY;
122  else {
123  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Must use mode 'standard', 'refmaxwell' or 'edge only'.");
124  }
125 
126  // If we're in edge only or standard modes, then the (2,2) hierarchy gets built without smoothers.
127  // Otherwise we use the user's smoothers (defaulting to Chebyshev if unspecified)
128  if(list.isSublist("maxwell1: 22list"))
129  precList22_ = list.sublist("maxwell1: 22list");
130  else if(list.isSublist("refmaxwell: 22list"))
131  precList22_ = list.sublist("refmaxwell: 22list");
132  if(mode_ == MODE_EDGE_ONLY || mode_ == MODE_STANDARD)
133  precList22_.set("smoother: pre or post","none");
134  else if(!precList22_.isType<std::string>("Preconditioner Type") &&
135  !precList22_.isType<std::string>("smoother: type") &&
136  !precList22_.isType<std::string>("smoother: pre type") &&
137  !precList22_.isType<std::string>("smoother: post type")) {
138  precList22_ = defaultSmootherList;
139  }
140  precList22_.set("verbosity",precList22_.get("verbosity",verbosity));
141 
142 
143  // For the (1,1) hierarchy we'll use Hiptmair (STANDARD) or Chevyshev (EDGE_ONLY / REFMAXWELL) if
144  // the user doesn't specify things
145  if(list.isSublist("maxwell1: 11list"))
146  precList11_ = list.sublist("maxwell1: 11list");
147  else if(list.isSublist("refmaxwell: 11list"))
148  precList11_ = list.sublist("refmaxwell: 11list");
149  if(!precList11_.isType<std::string>("Preconditioner Type") &&
150  !precList11_.isType<std::string>("smoother: type") &&
151  !precList11_.isType<std::string>("smoother: pre type") &&
152  !precList11_.isType<std::string>("smoother: post type")) {
153  if(mode_ == MODE_EDGE_ONLY || mode_ == MODE_REFMAXWELL) {
154  precList11_ = defaultSmootherList;
155  }
156  else {
157  precList11_.set("smoother: type", "HIPTMAIR");
158  precList11_.sublist("hiptmair: smoother type 1","CHEBYSHEV");
159  precList11_.sublist("hiptmair: smoother type 2","CHEBYSHEV");
160  precList11_.sublist("hiptmair: smoother list 1") = defaultSmootherList;
161  precList11_.sublist("hiptmair: smoother list 2") = defaultSmootherList;
162  }
163  }
164  precList11_.set("verbosity",precList11_.get("verbosity",verbosity));
165 
166  // Reuse support
167  if (enable_reuse_ &&
168  !precList11_.isType<std::string>("Preconditioner Type") &&
169  !precList11_.isParameter("reuse: type"))
170  precList11_.set("reuse: type", "full");
171  if (enable_reuse_ &&
172  !precList22_.isType<std::string>("Preconditioner Type") &&
173  !precList22_.isParameter("reuse: type"))
174  precList22_.set("reuse: type", "full");
175 
176 
177  // Are we using Kokkos?
178 #if !defined(HAVE_MUELU_KOKKOS_REFACTOR)
179  useKokkos_ = false;
180 #else
181 # ifdef HAVE_MUELU_SERIAL
182  if (typeid(Node).name() == typeid(Kokkos::Compat::KokkosSerialWrapperNode).name())
183  useKokkos_ = false;
184 # endif
185 # ifdef HAVE_MUELU_OPENMP
186  if (typeid(Node).name() == typeid(Kokkos::Compat::KokkosOpenMPWrapperNode).name())
187  useKokkos_ = true;
188 # endif
189 # ifdef HAVE_MUELU_CUDA
190  if (typeid(Node).name() == typeid(Kokkos::Compat::KokkosCudaWrapperNode).name())
191  useKokkos_ = true;
192 # endif
193  useKokkos_ = list.get("use kokkos refactor",useKokkos_);
194 #endif
195 
196  }
197 
198 
199  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
201 
202 #ifdef HAVE_MUELU_CUDA
203  if (parameterList_.get<bool>("maxwell1: cuda profile setup", false)) cudaProfilerStart();
204 #endif
205 
206  std::string timerLabel;
207  if (reuse)
208  timerLabel = "MueLu Maxwell1: compute (reuse)";
209  else
210  timerLabel = "MueLu Maxwell1: compute";
211  RCP<Teuchos::TimeMonitor> tmCompute = getTimer(timerLabel);
212 
214  // Remove explicit zeros from matrices
215  Maxwell_Utils<SC,LO,GO,NO>::removeExplicitZeros(parameterList_,D0_Matrix_,SM_Matrix_);
216 
217 
218  if (IsPrint(Statistics2)) {
219  RCP<ParameterList> params = rcp(new ParameterList());;
220  params->set("printLoadBalancingInfo", true);
221  params->set("printCommInfo", true);
222  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*SM_Matrix_, "SM_Matrix", params);
223  }
224 
226  // Detect Dirichlet boundary conditions
227  if (!reuse) {
228  magnitudeType rowSumTol = precList11_.get("aggregation: row sum drop tol",-1.0);
229  Maxwell_Utils<SC,LO,GO,NO>::detectBoundaryConditionsSM(SM_Matrix_,D0_Matrix_,rowSumTol,
231  useKokkos_,BCrowsKokkos_,BCcolsKokkos_,BCdomainKokkos_,
232 #endif
233  BCedges_,BCnodes_,BCrows_,BCcols_,BCdomain_,
234  allEdgesBoundary_,allNodesBoundary_);
235  if (IsPrint(Statistics2)) {
236  GetOStream(Statistics2) << "MueLu::Maxwell1::compute(): Detected " << BCedges_ << " BC rows and " << BCnodes_ << " BC columns." << std::endl;
237  }
238  }
239 
240  if (allEdgesBoundary_) {
241  // All edges have been detected as boundary edges.
242  // Do not attempt to construct sub-hierarchies, but just set up a single level preconditioner.
243  GetOStream(Warnings0) << "All edges are detected as boundary edges!" << std::endl;
244  mode_ = MODE_EDGE_ONLY;
245 
246  // Switch smoother off of Hiptmair
247  if(precList11_.get<std::string>("smoother: type") == "HIPTMAIR") {
248  precList11_.set("smoother: type",precList11_.get("hiptmair: smoother type 1","CHEBYSHEV"));
249  precList11_.sublist("smoother: sublist") = precList11_.sublist("hiptmair: smoother list 1");
250  }
251 
252  // Generate single level hierarchy for the edge
253  throw std::runtime_error("Maxwell1: Not yet supported");
254 
255  return;
256  }
257 
258  if (allNodesBoundary_) {
259  // All Nodes have been detected as boundary edges.
260  // Do not attempt to construct sub-hierarchies, but just set up a single level preconditioner.
261  GetOStream(Warnings0) << "All nodes are detected as boundary edges!" << std::endl;
262  mode_ = MODE_EDGE_ONLY;
263  // Switch smoother off of Hiptmair
264  if(precList11_.get<std::string>("smoother: type") == "HIPTMAIR") {
265  precList11_.set("smoother: type",precList11_.get("hiptmair: smoother type 1","CHEBYSHEV"));
266  precList11_.sublist("smoother: sublist") = precList11_.sublist("hiptmair: smoother list 1");
267  }
268 
269  // Generate single level hierarchy for the edge
270  throw std::runtime_error("Maxwell1: Not yet supported");
271 
272  return;
273  }
274 
275 
277  // Apply boundary conditions to D0 (if needed)
278  if(!reuse && applyBCsTo22_) {
279  GetOStream(Runtime0) << "Maxwell1::compute(): nuking BC nodes of D0" << std::endl;
280  D0_Matrix_->resumeFill();
281  Scalar replaceWith;
282  if (D0_Matrix_->getRowMap()->lib() == Xpetra::UseEpetra)
283  replaceWith= Teuchos::ScalarTraits<SC>::eps();
284  else
285  replaceWith = Teuchos::ScalarTraits<SC>::zero();
286 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
287  if (useKokkos_) {
288  Utilities_kokkos::ZeroDirichletCols(D0_Matrix_,BCcolsKokkos_,replaceWith);
289  } else {
290  Utilities::ZeroDirichletCols(D0_Matrix_,BCcols_,replaceWith);
291  }
292 #else
293  Utilities::ZeroDirichletCols(D0_Matrix_,BCcols_,replaceWith);
294 #endif
295  D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
296  }
297 
299  // Generate Kn and apply BCs (if needed)
300  if(Kn_Matrix_.is_null()) {
301  Kn_Matrix_ = generate_kn();
302  }
303 
304 
306  // Generate the (2,2) Hierarchy
307  Hierarchy22_ = MueLu::CreateXpetraPreconditioner(Kn_Matrix_, precList22_);
308 
309 
311  // Copy the relevant (2,2) data to the (1,1) hierarchy
312  Hierarchy11_ = rcp(new Hierarchy());
313  for(int i=0; i<Hierarchy22_->GetNumLevels(); i++) {
314  Hierarchy11_->AddNewLevel();
315  RCP<Level> NodeL = Hierarchy22_->GetLevel(i);
316  RCP<Level> EdgeL = Hierarchy11_->GetLevel(i);
317  EdgeL->Set("NodeMatrix",NodeL->Get<RCP<Matrix> >("A"));
318  if(i==0) {
319  EdgeL->Set("A", SM_Matrix_);
320  EdgeL->Set("D0", D0_Matrix_);
321  }
322  else {
323  EdgeL->Set("Pnodal",NodeL->Get<RCP<Matrix> >("P"));
324  }
325  }
326 
327 
329  // Generating the (1,1) Hierarchy
330  {
331  SM_Matrix_->setObjectLabel("A(1,1)");
332  precList11_.set("coarse: max size",1);
333  precList11_.set("max levels",Hierarchy22_->GetNumLevels());
334  // Hierarchy11_ = MueLu::CreateXpetraPreconditioner(SM_Matrix_, precList11_);
335  RCP<HierarchyManager<SC,LO,GO,NO> > mueLuFactory = rcp(new ParameterListInterpreter<SC,LO,GO,NO>(precList11_,SM_Matrix_->getDomainMap()->getComm()));
336  Hierarchy11_->setlib(SM_Matrix_->getDomainMap()->lib());
337  Hierarchy11_->SetProcRankVerbose(SM_Matrix_->getDomainMap()->getComm()->getRank());
338  mueLuFactory->SetupHierarchy(*Hierarchy11_);
339 
340  if(mode_ == MODE_REFMAXWELL) {
341  if(Hierarchy11_->GetNumLevels() > 1) {
342  RCP<Level> EdgeL = Hierarchy11_->GetLevel(1);
343  P11_ = EdgeL->Get<RCP<Matrix> >("P");
344  }
345  }
346  }
347 
349  // Allocate temporary MultiVectors for solve (only needed for RefMaxwell style)
350  allocateMemory(1);
351 
352  }
353 
354 
355  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
356  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::generate_kn() const {
357  // NOTE: This does not nicely support reuse, but the relevant coed can be copied from
358  // RefMaxwell when we decide we want to do this.
359 
360  // NOTE: Boundary conditions OAZ are handled via the "rap: fix zero diagonals threshold"
361  RCP<Teuchos::TimeMonitor> tm = getTimer("MueLu Maxwell1: Build Kn");
362 
363  Level fineLevel, coarseLevel;
364  fineLevel.SetFactoryManager(null);
365  coarseLevel.SetFactoryManager(null);
366  coarseLevel.SetPreviousLevel(rcpFromRef(fineLevel));
367  fineLevel.SetLevelID(0);
368  coarseLevel.SetLevelID(1);
369  fineLevel.Set("A",SM_Matrix_);
370  coarseLevel.Set("P",D0_Matrix_);
371  //coarseLevel.Set("Coordinates",Coords_);
372 
373  coarseLevel.setlib(SM_Matrix_->getDomainMap()->lib());
374  fineLevel.setlib(SM_Matrix_->getDomainMap()->lib());
375  coarseLevel.setObjectLabel("Maxwell1 (2,2)");
376  fineLevel.setObjectLabel("Maxwell1 (2,2)");
377 
378  RCP<RAPFactory> rapFact = rcp(new RAPFactory());
379  ParameterList rapList = *(rapFact->GetValidParameterList());
380  rapList.set("transpose: use implicit", true);
381  rapList.set("rap: fix zero diagonals", parameterList_.get<bool>("rap: fix zero diagonals", true));
382  rapList.set("rap: fix zero diagonals threshold", parameterList_.get<double>("rap: fix zero diagonals threshold", Teuchos::ScalarTraits<double>::eps()));
383  rapList.set("rap: triple product", parameterList_.get<bool>("rap: triple product", false));
384  rapFact->SetParameterList(rapList);
385  coarseLevel.Request("A", rapFact.get());
386  if (enable_reuse_) {
387  coarseLevel.Request("AP reuse data", rapFact.get());
388  coarseLevel.Request("RAP reuse data", rapFact.get());
389  }
390 
391  RCP<Matrix> Kn_Matrix = coarseLevel.Get< RCP<Matrix> >("A", rapFact.get());
392  Kn_Matrix->setObjectLabel("A(2,2)");
393 
394  return Kn_Matrix;
395  }
396 
397 
398 
399  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
401  if(mode_ == MODE_REFMAXWELL) {
402  RCP<Teuchos::TimeMonitor> tmAlloc = getTimer("MueLu Maxwell1: Allocate MVs");
403 
404  residualFine_ = MultiVectorFactory::Build(SM_Matrix_->getRangeMap(), numVectors);
405 
406  if(!Hierarchy11_.is_null() && Hierarchy11_->GetNumLevels() > 1) {
407  RCP<Level> EdgeL = Hierarchy11_->GetLevel(1);
408  RCP<Matrix> A = EdgeL->Get<RCP<Matrix> >("A");
409  residual11c_ = MultiVectorFactory::Build(A->getRangeMap(), numVectors);
410  update11c_ = MultiVectorFactory::Build(A->getDomainMap(), numVectors);
411  }
412 
413  if(!Hierarchy22_.is_null()) {
414  residual22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
415  update22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
416  }
417 
418  }
419 
420  }
421 
422 
423  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
424  void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dump(const Matrix& A, std::string name) const {
425  if (dump_matrices_) {
426  GetOStream(Runtime0) << "Dumping to " << name << std::endl;
427  Xpetra::IO<SC, LO, GO, NO>::Write(name, A);
428  }
429  }
430 
431 
432  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
433  void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dump(const MultiVector& X, std::string name) const {
434  if (dump_matrices_) {
435  GetOStream(Runtime0) << "Dumping to " << name << std::endl;
436  Xpetra::IO<SC, LO, GO, NO>::Write(name, X);
437  }
438  }
439 
440 
441  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
443  if (dump_matrices_) {
444  GetOStream(Runtime0) << "Dumping to " << name << std::endl;
445  Xpetra::IO<coordinateType, LO, GO, NO>::Write(name, X);
446  }
447  }
448 
449 
450  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
451  void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dump(const Teuchos::ArrayRCP<bool>& v, std::string name) const {
452  if (dump_matrices_) {
453  GetOStream(Runtime0) << "Dumping to " << name << std::endl;
454  std::ofstream out(name);
455  for (size_t i = 0; i < Teuchos::as<size_t>(v.size()); i++)
456  out << v[i] << "\n";
457  }
458  }
459 
460 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
461  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
463  if (dump_matrices_) {
464  GetOStream(Runtime0) << "Dumping to " << name << std::endl;
465  std::ofstream out(name);
466  auto vH = Kokkos::create_mirror_view (v);
467  Kokkos::deep_copy(vH , v);
468  for (size_t i = 0; i < vH.size(); i++)
469  out << vH[i] << "\n";
470  }
471  }
472 #endif
473 
474  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
475  Teuchos::RCP<Teuchos::TimeMonitor> Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getTimer(std::string name, RCP<const Teuchos::Comm<int> > comm) const {
476  if (IsPrint(Timings)) {
477  if (!syncTimers_)
478  return Teuchos::rcp(new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name)));
479  else {
480  if (comm.is_null())
481  return Teuchos::rcp(new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), SM_Matrix_->getRowMap()->getComm().ptr()));
482  else
483  return Teuchos::rcp(new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), comm.ptr()));
484  }
485  } else
486  return Teuchos::null;
487  }
488 
489 
490 
491 
492  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
493  void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::resetMatrix(RCP<Matrix> SM_Matrix_new, bool ComputePrec) {
494  bool reuse = !SM_Matrix_.is_null();
495  SM_Matrix_ = SM_Matrix_new;
496  dump(*SM_Matrix_, "SM.m");
497  if (ComputePrec)
498  compute(reuse);
499  }
500 
501 
502  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
503  void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::applyInverseRefMaxwellAdditive(const MultiVector& RHS, MultiVector& X) const {
504  // make sure that we have enough temporary memory
505  const SC one = Teuchos::ScalarTraits<SC>::one();
506  if (!allEdgesBoundary_ && X.getNumVectors() != residualFine_->getNumVectors())
507  allocateMemory(X.getNumVectors());
508 
509  TEUCHOS_TEST_FOR_EXCEPTION(Hierarchy11_.is_null() || Hierarchy11_->GetNumLevels() == 0, Exceptions::RuntimeError, "(1,1) Hiearchy is null.");
510 
511  // 1) Run fine pre-smoother using Hierarchy11
512  RCP<Level> Fine = Hierarchy11_->GetLevel(0);
513  if (Fine->IsAvailable("PreSmoother")) {
514  RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: PreSmoother");
515  RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
516  preSmoo->Apply(X, RHS, true);
517  }
518 
519  // 2) Compute residual
520  {
521  RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: residual calculation");
522  Utilities::Residual(*SM_Matrix_, X, RHS,*residualFine_);
523  }
524 
525  // 3a) Restrict residual to (1,1) Hierarchy's level 1 and execute (1,1) hierarchy (use startLevel and InitialGuessIsZero)
526  if(!P11_.is_null()) {
527  RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: (1,1) correction");
528  P11_->apply(*residualFine_,*residual11c_,Teuchos::TRANS);
529  Hierarchy11_->Iterate(*residual11c_,*update11c_,true,1);
530  }
531 
532  // 3b) Restrict residual to (2,2) Hierarchy's level 0 and execute (2,2) hierarchy (use InitialGuessIsZero)
533  if (!allNodesBoundary_) {
534  RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: (2,2) correction");
535  D0_Matrix_->apply(*residualFine_,*residual22_,Teuchos::TRANS);
536  Hierarchy22_->Iterate(*residual22_,*update22_,true,0);
537  }
538 
539  // 4) Prolong both updates back into X-vector (Need to do both the P11 null and not null cases
540  {
541  RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: Orolongation");
542  if(!P11_.is_null())
543  P11_->apply(*update11c_,X,Teuchos::NO_TRANS,one,one);
544  if (!allNodesBoundary_)
545  D0_Matrix_->apply(*update22_,X,Teuchos::NO_TRANS,one,one);
546  }
547 
548  // 5) Run fine post-smoother using Hierarchy11
549  if (Fine->IsAvailable("PostSmoother")) {
550  RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: PostSmoother");
551  RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
552  postSmoo->Apply(X, RHS, false);
553  }
554 
555 
556  }
557 
558  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
559  void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::applyInverseStandard(const MultiVector& RHS, MultiVector& X) const {
560  Hierarchy11_->Iterate(RHS,X);
561  }
562 
563  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
564  void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::apply (const MultiVector& RHS, MultiVector& X,
565  Teuchos::ETransp /* mode */,
566  Scalar /* alpha */,
567  Scalar /* beta */) const {
568  RCP<Teuchos::TimeMonitor> tm = getTimer("MueLu Maxwell1: solve");
569  if(mode_ == MODE_STANDARD || mode_ == MODE_EDGE_ONLY)
570  applyInverseStandard(RHS,X);
571  else if(mode_ == MODE_REFMAXWELL)
572  applyInverseRefMaxwellAdditive(RHS,X);
573  else
574  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Must use mode 'standard', 'refmaxwell' or 'edge only'.");
575  }
576 
577 
578 
579  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
581  return false;
582  }
583 
584 
585  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
587  initialize(const Teuchos::RCP<Matrix> & D0_Matrix,
588  const Teuchos::RCP<Matrix> & Kn_Matrix,
589  const Teuchos::RCP<MultiVector> & Nullspace,
590  const Teuchos::RCP<RealValuedMultiVector> & Coords,
591  Teuchos::ParameterList& List)
592  {
593  // some pre-conditions
594  TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
595 
596 #ifdef HAVE_MUELU_DEBUG
597  if(!Kn_Matrix.is_null()) {
598  TEUCHOS_ASSERT(Kn_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getDomainMap()));
599  TEUCHOS_ASSERT(Kn_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getDomainMap()));
600  }
601 
602 
603  TEUCHOS_ASSERT(D0_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getRowMap()));
604 #endif
605 
606  Hierarchy11_ = Teuchos::null;
607  Hierarchy22_ = Teuchos::null;
608  mode_ = MODE_STANDARD;
609 
610  // Default settings
611  useKokkos_=false;
612  allEdgesBoundary_=false;
613  allNodesBoundary_=false;
614  dump_matrices_ = false;
615  enable_reuse_=false;
616  syncTimers_=false;
617  applyBCsTo22_ = true;
618 
619 
620 
621  // set parameters
622  setParameters(List);
623 
624  if (D0_Matrix->getRowMap()->lib() == Xpetra::UseTpetra) {
625  // We will remove boundary conditions from D0, and potentially change maps, so we copy the input.
626  // Fortunately, D0 is quite sparse.
627  // We cannot use the Tpetra copy constructor, since it does not copy the graph.
628 
629  RCP<Matrix> D0copy = MatrixFactory::Build(D0_Matrix->getRowMap(), D0_Matrix->getColMap(), 0);
630  RCP<CrsMatrix> D0copyCrs = rcp_dynamic_cast<CrsMatrixWrap>(D0copy,true)->getCrsMatrix();
631  ArrayRCP<const size_t> D0rowptr_RCP;
632  ArrayRCP<const LO> D0colind_RCP;
633  ArrayRCP<const SC> D0vals_RCP;
634  rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,true)->getCrsMatrix()->getAllValues(D0rowptr_RCP,
635  D0colind_RCP,
636  D0vals_RCP);
637 
638  ArrayRCP<size_t> D0copyrowptr_RCP;
639  ArrayRCP<LO> D0copycolind_RCP;
640  ArrayRCP<SC> D0copyvals_RCP;
641  D0copyCrs->allocateAllValues(D0vals_RCP.size(),D0copyrowptr_RCP,D0copycolind_RCP,D0copyvals_RCP);
642  D0copyrowptr_RCP.deepCopy(D0rowptr_RCP());
643  D0copycolind_RCP.deepCopy(D0colind_RCP());
644  D0copyvals_RCP.deepCopy(D0vals_RCP());
645  D0copyCrs->setAllValues(D0copyrowptr_RCP,
646  D0copycolind_RCP,
647  D0copyvals_RCP);
648  D0copyCrs->expertStaticFillComplete(D0_Matrix->getDomainMap(), D0_Matrix->getRangeMap(),
649  rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,true)->getCrsMatrix()->getCrsGraph()->getImporter(),
650  rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,true)->getCrsMatrix()->getCrsGraph()->getExporter());
651  D0_Matrix_ = D0copy;
652  } else
653  D0_Matrix_ = MatrixFactory::BuildCopy(D0_Matrix);
654 
655 
656  Kn_Matrix_ = Kn_Matrix;
657  Coords_ = Coords;
658  Nullspace_ = Nullspace;
659 
660  if(!Kn_Matrix_.is_null()) dump(*Kn_Matrix_, "Kn.m");
661  if (!Nullspace_.is_null()) dump(*Nullspace_, "nullspace.m");
662  if (!Coords_.is_null()) dumpCoords(*Coords_, "coords.m");
663 
664  }
665 
666 
667 
668  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
670  describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel /* verbLevel */) const {
671 
672  std::ostringstream oss;
673 
674  RCP<const Teuchos::Comm<int> > comm = SM_Matrix_->getDomainMap()->getComm();
675 
676 #ifdef HAVE_MPI
677  int root;
678  if (!Kn_Matrix_.is_null())
679  root = comm->getRank();
680  else
681  root = -1;
682 
683  int actualRoot;
684  reduceAll(*comm, Teuchos::REDUCE_MAX, root, Teuchos::ptr(&actualRoot));
685  root = actualRoot;
686 #endif
687 
688 
689  oss << "\n--------------------------------------------------------------------------------\n" <<
690  "--- Maxwell1 Summary ---\n"
691  "--------------------------------------------------------------------------------" << std::endl;
692  oss << std::endl;
693 
694  GlobalOrdinal numRows;
695  GlobalOrdinal nnz;
696 
697  SM_Matrix_->getRowMap()->getComm()->barrier();
698 
699  numRows = SM_Matrix_->getGlobalNumRows();
700  nnz = SM_Matrix_->getGlobalNumEntries();
701 
702  Xpetra::global_size_t tt = numRows;
703  int rowspacer = 3; while (tt != 0) { tt /= 10; rowspacer++; }
704  tt = nnz;
705  int nnzspacer = 2; while (tt != 0) { tt /= 10; nnzspacer++; }
706 
707  oss << "block " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz " << std::setw(9) << " nnz/row" << std::endl;
708  oss << "(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
709 
710  if (!Kn_Matrix_.is_null()) {
711  // ToDo: make sure that this is printed correctly
712  numRows = Kn_Matrix_->getGlobalNumRows();
713  nnz = Kn_Matrix_->getGlobalNumEntries();
714 
715  oss << "(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
716  }
717 
718  oss << std::endl;
719 
720 
721  std::string outstr = oss.str();
722 
723 #ifdef HAVE_MPI
724  RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
725  MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
726 
727  int strLength = outstr.size();
728  MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
729  if (comm->getRank() != root)
730  outstr.resize(strLength);
731  MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
732 #endif
733 
734  out << outstr;
735 
736  if (!Hierarchy11_.is_null())
737  Hierarchy11_->describe(out, GetVerbLevel());
738 
739  if (!Hierarchy22_.is_null())
740  Hierarchy22_->describe(out, GetVerbLevel());
741 
742  if (IsPrint(Statistics2)) {
743  // Print the grid of processors
744  std::ostringstream oss2;
745 
746  oss2 << "Sub-solver distribution over ranks" << std::endl;
747  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;
748 
749  int numProcs = comm->getSize();
750 #ifdef HAVE_MPI
751  RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
752 
753  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
754 #endif
755 
756  char status = 0;
757  if (!Kn_Matrix_.is_null())
758  status += 1;
759  std::vector<char> states(numProcs, 0);
760 #ifdef HAVE_MPI
761  MPI_Gather(&status, 1, MPI_CHAR, &states[0], 1, MPI_CHAR, 0, *rawMpiComm);
762 #else
763  states.push_back(status);
764 #endif
765 
766  int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
767  for (int proc = 0; proc < numProcs; proc += rowWidth) {
768  for (int j = 0; j < rowWidth; j++)
769  if (proc + j < numProcs)
770  if (states[proc+j] == 0)
771  oss2 << ".";
772  else if (states[proc+j] == 1)
773  oss2 << "1";
774  else if (states[proc+j] == 2)
775  oss2 << "2";
776  else
777  oss2 << "B";
778  else
779  oss2 << " ";
780 
781  oss2 << " " << proc << ":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
782  }
783  oss2 << std::endl;
784  GetOStream(Statistics2) << oss2.str();
785  }
786 
787 
788  }
789 
790 
791 
792 } // namespace
793 
794 #define MUELU_MAXWELL1_SHORT
795 #endif //ifdef MUELU_MAXWELL1_DEF_HPP
Important warning messages (one line)
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
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.
static void removeExplicitZeros(Teuchos::ParameterList &parameterList, RCP< Matrix > &D0_Matrix, RCP< Matrix > &SM_Matrix, RCP< Matrix > &M1_Matrix, RCP< Matrix > &Ms_Matrix)
Remove explicit zeros.
Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
void applyInverseStandard(const MultiVector &RHS, MultiVector &X) const
apply standard Maxwell1 cycle
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
void applyInverseRefMaxwellAdditive(const MultiVector &RHS, MultiVector &X) const
apply RefMaxwell additive 2x2 style cycle
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
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)
One-liner description of what is happening.
Namespace for MueLu classes and methods.
void SetPreviousLevel(const RCP< Level > &previousLevel)
Definition: MueLu_Level.cpp:85
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
Definition: MueLu_Level.cpp:92
MueLu::DefaultNode Node
Print even more statistics.
void setlib(Xpetra::UnderlyingLib lib2)
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
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)
void setParameters(Teuchos::ParameterList &list)
Set parameters.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Teuchos::RCP< Teuchos::TimeMonitor > getTimer(std::string name, RCP< const Teuchos::Comm< int > > comm=Teuchos::null) const
get a (synced) timer
static void detectBoundaryConditionsSM(RCP< Matrix > &SM_Matrix, RCP< Matrix > &D0_Matrix, magnitudeType rowSumTol, 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 std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Print all timing information.
void dump(const Matrix &A, std::string name) const
dump out matrix
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
void SetLevelID(int levelID)
Set level number.
Definition: MueLu_Level.cpp:78
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...
void compute(bool reuse=false)
Setup the preconditioner.
Teuchos::RCP< Matrix > generate_kn() const
Generates the Kn matrix.
Xpetra::MultiVector< coordinateType, LO, GO, NO > RealValuedMultiVector
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.
static void ZeroDirichletCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
Factory for building coarse matrices.
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
void dumpCoords(const RealValuedMultiVector &X, std::string name) const
dump out real-valued multivector
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
void allocateMemory(int numVectors) const
allocate multivectors for solve