MueLu Version of the Day
Loading...
Searching...
No Matches
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_CrsMatrixUtils.hpp"
54#include "Xpetra_MatrixUtils.hpp"
55
57#include "MueLu_Maxwell_Utils.hpp"
58
59#include "MueLu_ReitzingerPFactory.hpp"
60#include "MueLu_SaPFactory.hpp"
61#include "MueLu_AggregationExportFactory.hpp"
62#include "MueLu_Utilities.hpp"
63#include "MueLu_Level.hpp"
64#include "MueLu_Hierarchy.hpp"
65#include "MueLu_RAPFactory.hpp"
66#include "MueLu_Monitor.hpp"
67#include "MueLu_PerfUtils.hpp"
68#include "MueLu_ParameterListInterpreter.hpp"
70#include <MueLu_HierarchyUtils.hpp>
71# include "MueLu_Utilities_kokkos.hpp"
75#include <MueLu_RefMaxwellSmoother.hpp>
76
77#ifdef HAVE_MUELU_CUDA
78#include "cuda_profiler_api.h"
79#endif
80
81
82namespace MueLu {
83
84
85 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getDomainMap() const {
87 return SM_Matrix_->getDomainMap();
88 }
89
90
91 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
92 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getRangeMap() const {
93 return SM_Matrix_->getRangeMap();
94 }
95
96
97 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
99
100 if (list.isType<std::string>("parameterlist: syntax") && list.get<std::string>("parameterlist: syntax") == "ml") {
101 list.remove("parameterlist: syntax");
102 Teuchos::ParameterList newList;
103
104 // interpret ML list
105 newList.sublist("maxwell1: 22list") = *Teuchos::getParametersFromXmlString(MueLu::ML2MueLuParameterTranslator::translate(list,"Maxwell"));
106
107
108 // Hardwiring options to ensure ML compatibility
109 newList.sublist("maxwell1: 22list").set("use kokkos refactor", false);
110 newList.sublist("maxwell1: 22list").set("tentative: constant column sums", false);
111 newList.sublist("maxwell1: 22list").set("tentative: calculate qr", false);
112
113 newList.sublist("maxwell1: 11list").set("use kokkos refactor", false);
114 newList.sublist("maxwell1: 11list").set("tentative: constant column sums", false);
115 newList.sublist("maxwell1: 11list").set("tentative: calculate qr", false);
116
117 if(list.isParameter("aggregation: damping factor") && list.get<double>("aggregation: damping factor") == 0.0)
118 newList.sublist("maxwell1: 11list").set("multigrid algorithm", "unsmoothed reitzinger");
119 else
120 newList.sublist("maxwell1: 11list").set("multigrid algorithm", "smoothed reitzinger");
121 newList.sublist("maxwell1: 11list").set("aggregation: type", "uncoupled");
122
123 newList.sublist("maxwell1: 22list").set("multigrid algorithm", "unsmoothed");
124 newList.sublist("maxwell1: 22list").set("aggregation: type", "uncoupled");
125
126 if (newList.sublist("maxwell1: 22list").isType<std::string>("verbosity"))
127 newList.set("verbosity", newList.sublist("maxwell1: 22list").get<std::string>("verbosity"));
128
129 // Move coarse solver and smoother stuff to 11list
130 std::vector<std::string> convert = {"coarse:", "smoother:", "smoother: pre", "smoother: post"};
131 for (auto it = convert.begin(); it != convert.end(); ++it) {
132 if (newList.sublist("maxwell1: 22list").isType<std::string>(*it + " type")) {
133 newList.sublist("maxwell1: 11list").set(*it+" type", newList.sublist("maxwell1: 22list").get<std::string>(*it+" type"));
134 newList.sublist("maxwell1: 22list").remove(*it+" type");
135 }
136 if (newList.sublist("maxwell1: 22list").isSublist(*it+" params")) {
137 newList.sublist("maxwell1: 11list").set(*it+" params", newList.sublist("maxwell1: 22list").sublist(*it+" params"));
138 newList.sublist("maxwell1: 22list").remove(*it+" params");
139 }
140 }
141
142 newList.sublist("maxwell1: 22list").set("smoother: type", "none");
143 newList.sublist("maxwell1: 22list").set("coarse: type", "none");
144
145 list = newList;
146 }
147 std::string mode_string = list.get("maxwell1: mode", MasterList::getDefault<std::string>("maxwell1: mode"));
148 applyBCsTo22_ = list.get("maxwell1: apply BCs to 22", false);
149 dump_matrices_ = list.get("maxwell1: dump matrices", MasterList::getDefault<bool>("maxwell1: dump matrices"));
150
151 // Default smoother. We'll copy this around.
152 Teuchos::ParameterList defaultSmootherList;
153 defaultSmootherList.set("smoother: type", "CHEBYSHEV");
154 defaultSmootherList.sublist("smoother: params").set("chebyshev: degree",2);
155 defaultSmootherList.sublist("smoother: params").set("chebyshev: ratio eigenvalue",7.0);
156 defaultSmootherList.sublist("smoother: params").set("chebyshev: eigenvalue max iterations",30);
157
158 // Make sure verbosity gets passed to the sublists
159 std::string verbosity = list.get("verbosity","high");
161
162 // Check the validity of the run mode
163 if(mode_ != MODE_GMHD_STANDARD) {
164 if(mode_string == "standard") mode_ = MODE_STANDARD;
165 else if(mode_string == "refmaxwell") mode_ = MODE_REFMAXWELL;
166 else if(mode_string == "edge only") mode_ = MODE_EDGE_ONLY;
167 else {
168 TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Must use mode 'standard', 'refmaxwell', 'edge only', or use the GMHD constructor.");
169 }
170 }
171
172 // If we're in edge only or standard modes, then the (2,2) hierarchy gets built without smoothers.
173 // Otherwise we use the user's smoothers (defaulting to Chebyshev if unspecified)
174 if(list.isSublist("maxwell1: 22list"))
175 precList22_ = list.sublist("maxwell1: 22list");
176 else if(list.isSublist("refmaxwell: 22list"))
177 precList22_ = list.sublist("refmaxwell: 22list");
178 if(mode_ == MODE_EDGE_ONLY || mode_ == MODE_STANDARD || mode_ == MODE_GMHD_STANDARD)
179 precList22_.set("smoother: pre or post","none");
180 else if(!precList22_.isType<std::string>("Preconditioner Type") &&
181 !precList22_.isType<std::string>("smoother: type") &&
182 !precList22_.isType<std::string>("smoother: pre type") &&
183 !precList22_.isType<std::string>("smoother: post type")) {
184 precList22_ = defaultSmootherList;
185 }
186 precList22_.set("verbosity",precList22_.get("verbosity",verbosity));
187
188
189
190 // For the (1,1) hierarchy we'll use Hiptmair (STANDARD) or Chebyshev (EDGE_ONLY / REFMAXWELL) if
191 // the user doesn't specify things
192 if(list.isSublist("maxwell1: 11list"))
193 precList11_ = list.sublist("maxwell1: 11list");
194 else if(list.isSublist("refmaxwell: 11list"))
195 precList11_ = list.sublist("refmaxwell: 11list");
196
197 if(mode_ == MODE_GMHD_STANDARD) {
198 precList11_.set("smoother: pre or post","none");
199 precList11_.set("smoother: type", "none");
200 }
201 if(!precList11_.isType<std::string>("Preconditioner Type") &&
202 !precList11_.isType<std::string>("smoother: type") &&
203 !precList11_.isType<std::string>("smoother: pre type") &&
204 !precList11_.isType<std::string>("smoother: post type")) {
205 if(mode_ == MODE_EDGE_ONLY || mode_ == MODE_REFMAXWELL) {
206 precList11_ = defaultSmootherList;
207 }
208 if (mode_ == MODE_STANDARD) {
209 precList11_.set("smoother: type", "HIPTMAIR");
210 precList11_.sublist("hiptmair: smoother type 1","CHEBYSHEV");
211 precList11_.sublist("hiptmair: smoother type 2","CHEBYSHEV");
212 precList11_.sublist("hiptmair: smoother list 1") = defaultSmootherList;
213 precList11_.sublist("hiptmair: smoother list 2") = defaultSmootherList;
214 }
215 }
216 precList11_.set("verbosity",precList11_.get("verbosity",verbosity));
217
218 // Reuse support
219 if (enable_reuse_ &&
220 !precList11_.isType<std::string>("Preconditioner Type") &&
221 !precList11_.isParameter("reuse: type"))
222 precList11_.set("reuse: type", "full");
223 if (enable_reuse_ &&
224 !precList22_.isType<std::string>("Preconditioner Type") &&
225 !precList22_.isParameter("reuse: type"))
226 precList22_.set("reuse: type", "full");
227
228
229 // Are we using Kokkos?
230# ifdef HAVE_MUELU_SERIAL
231 if (typeid(Node).name() == typeid(Kokkos::Compat::KokkosSerialWrapperNode).name())
232 useKokkos_ = false;
233# endif
234# ifdef HAVE_MUELU_OPENMP
235 if (typeid(Node).name() == typeid(Kokkos::Compat::KokkosOpenMPWrapperNode).name())
236 useKokkos_ = true;
237# endif
238# ifdef HAVE_MUELU_CUDA
239 if (typeid(Node).name() == typeid(Kokkos::Compat::KokkosCudaWrapperNode).name())
240 useKokkos_ = true;
241# endif
242 useKokkos_ = list.get("use kokkos refactor",useKokkos_);
243
244 parameterList_ = list;
245
246 }
247
248 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
250 Teuchos::ParameterList precListGmhd;
251
252 MueLu::HierarchyUtils<SC,LO,GO,NO>::CopyBetweenHierarchies(*Hierarchy11_,*HierarchyGmhd_, "P", "Psubblock", "RCP<Matrix>");
253
254 HierarchyGmhd_->GetLevel(0)->Set("A", GmhdA_Matrix_);
255 GmhdA_Matrix_->setObjectLabel("GmhdA");
256
257 TEUCHOS_TEST_FOR_EXCEPTION( !List.isSublist("maxwell1: Gmhdlist"), Exceptions::RuntimeError, "Must provide maxwell1: Gmhdlist for GMHD setup");
258 precListGmhd = List.sublist("maxwell1: Gmhdlist");
259 precListGmhd.set("coarse: max size",1);
260 precListGmhd.set("max levels",HierarchyGmhd_->GetNumLevels());
261 RCP<MueLu::HierarchyManager<SC,LO,GO,NO> > mueLuFactory = rcp(new MueLu::ParameterListInterpreter<SC,LO,GO,NO>(precListGmhd,GmhdA_Matrix_->getDomainMap()->getComm()));
262 HierarchyGmhd_->setlib(GmhdA_Matrix_->getDomainMap()->lib());
263 HierarchyGmhd_->SetProcRankVerbose(GmhdA_Matrix_->getDomainMap()->getComm()->getRank());
264 mueLuFactory->SetupHierarchy(*HierarchyGmhd_);
265 }
266
267 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
269 /* Algorithm overview for Maxwell1 construction:
270
271 1) Create a nodal auxillary hierarchy based on (a) the user's nodal matrix or (b) a matrix constructed
272 by D0^T A D0 if the user doesn't provide a nodal matrix. We call this matrix "NodeAggMatrix."
273
274 2) If the user provided a node matrix, we use the prolongators from the auxillary nodal hierarchy
275 to generate matrices for smoothers on all levels. We call this "NodeMatrix." Otherwise NodeMatrix = NodeAggMatrix
276
277 3) We stick all of the nodal P matrices and NodeMatrix objects on the final (1,1) hierarchy and then use the
278 ReitzingerPFactory to generate the edge P and A matrices.
279 */
280
281
282#ifdef HAVE_MUELU_CUDA
283 if (parameterList_.get<bool>("maxwell1: cuda profile setup", false)) cudaProfilerStart();
284#endif
285
286 std::string timerLabel;
287 if (reuse)
288 timerLabel = "MueLu Maxwell1: compute (reuse)";
289 else
290 timerLabel = "MueLu Maxwell1: compute";
291 RCP<Teuchos::TimeMonitor> tmCompute = getTimer(timerLabel);
292
294 // Generate Kn and apply BCs (if needed)
295 bool have_generated_Kn = false;
296 if(Kn_Matrix_.is_null()) {
297 GetOStream(Runtime0) << "Maxwell1::compute(): Kn not provided. Generating." << std::endl;
298 Kn_Matrix_ = generate_kn();
299 have_generated_Kn = true;
300 }
301
302 if (parameterList_.get<bool>("rap: fix zero diagonals", true)) {
303 magnitudeType threshold;
304 if (parameterList_.isType<magnitudeType>("rap: fix zero diagonals threshold"))
305 threshold = parameterList_.get<magnitudeType>("rap: fix zero diagonals threshold",
306 Teuchos::ScalarTraits<double>::eps());
307 else
308 threshold = Teuchos::as<magnitudeType>(parameterList_.get<double>("rap: fix zero diagonals threshold",
309 Teuchos::ScalarTraits<double>::eps()));
310 Scalar replacement = Teuchos::as<Scalar>(parameterList_.get<double>("rap: fix zero diagonals replacement",
311 MasterList::getDefault<double>("rap: fix zero diagonals replacement")));
312 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Kn_Matrix_, true, GetOStream(Warnings1), threshold, replacement);
313 }
314
315
317 // Generate the (2,2) Hierarchy
318 Kn_Matrix_->setObjectLabel("Maxwell1 (2,2)");
319
320 /* Critical ParameterList changes */
321 if (!Coords_.is_null())
322 precList22_.sublist("user data").set("Coordinates",Coords_);
323
324 /* Repartitioning *must* be in sync between hierarchies, but the
325 only thing we need to watch here is the subcomms, since ReitzingerPFactory
326 won't look at all the other stuff */
327 if(precList22_.isParameter("repartition: enable")) {
328 bool repartition = precList22_.get<bool>("repartition: enable");
329 precList11_.set("repartition: enable",repartition);
330
331 // If we're repartitioning (2,2), we need to rebalance for (1,1) to do the right thing
332 if(repartition)
333 precList22_.set("repartition: rebalance P and R",true);
334
335 if (precList22_.isParameter("repartition: use subcommunicators")) {
336 precList11_.set("repartition: use subcommunicators", precList22_.get<bool>("repartition: use subcommunicators"));
337
338 // We do not want (1,1) and (2,2) blocks being repartitioned seperately, so we specify the map that
339 // is going to be used (this is generated in ReitzingerPFactory)
340 if(precList11_.get<bool>("repartition: use subcommunicators")==true)
341 precList11_.set("repartition: use subcommunicators in place",true);
342 }
343 else {
344 // We'll have Maxwell1 default to using subcommunicators if you don't specify
345 precList11_.set("repartition: use subcommunicators", true);
346 precList22_.set("repartition: use subcommunicators", true);
347
348 // We do not want (1,1) and (2,2) blocks being repartitioned seperately, so we specify the map that
349 // is going to be used (this is generated in ReitzingerPFactory)
350 precList11_.set("repartition: use subcommunicators in place",true);
351 }
352
353 }
354 else
355 precList11_.remove("repartition: enable", false);
356
357
359 // Remove explicit zeros from matrices
360 /*
361 Maxwell_Utils<SC,LO,GO,NO>::removeExplicitZeros(parameterList_,D0_Matrix_,SM_Matrix_);
362
363
364 if (IsPrint(Statistics2)) {
365 RCP<ParameterList> params = rcp(new ParameterList());;
366 params->set("printLoadBalancingInfo", true);
367 params->set("printCommInfo", true);
368 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*SM_Matrix_, "SM_Matrix", params);
369 }
370 */
371
372
374 // Detect Dirichlet boundary conditions
375 if (!reuse) {
376 magnitudeType rowSumTol = precList11_.get("aggregation: row sum drop tol",-1.0);
377 Maxwell_Utils<SC,LO,GO,NO>::detectBoundaryConditionsSM(SM_Matrix_,D0_Matrix_,rowSumTol,
378 useKokkos_,BCrowsKokkos_,BCcolsKokkos_,BCdomainKokkos_,
379 BCedges_,BCnodes_,BCrows_,BCcols_,BCdomain_,
380 allEdgesBoundary_,allNodesBoundary_);
381 if (IsPrint(Statistics2)) {
382 GetOStream(Statistics2) << "MueLu::Maxwell1::compute(): Detected " << BCedges_ << " BC rows and " << BCnodes_ << " BC columns." << std::endl;
383 }
384 }
385
386 if (allEdgesBoundary_) {
387 // All edges have been detected as boundary edges.
388 // Do not attempt to construct sub-hierarchies, but just set up a single level preconditioner.
389 GetOStream(Warnings0) << "All edges are detected as boundary edges!" << std::endl;
390 mode_ = MODE_EDGE_ONLY;
391
392 // Generate single level hierarchy for the edge
393 precList22_.set("max levels", 1);
394 }
395
396
397 if (allNodesBoundary_) {
398 // All Nodes have been detected as boundary nodes.
399 // Do not attempt to construct sub-hierarchies, but just set up a single level preconditioner.
400 GetOStream(Warnings0) << "All nodes are detected as boundary nodes!" << std::endl;
401 mode_ = MODE_EDGE_ONLY;
402
403 // Generate single level hierarchy for the edge
404 precList22_.set("max levels", 1);
405 }
406
408 // Build (2,2) hierarchy
409 Hierarchy22_ = MueLu::CreateXpetraPreconditioner(Kn_Matrix_, precList22_);
410
411
413 // Apply boundary conditions to D0 (if needed)
414 if(!reuse) {
415 D0_Matrix_->resumeFill();
416 Scalar replaceWith;
417 if (D0_Matrix_->getRowMap()->lib() == Xpetra::UseEpetra)
418 replaceWith= Teuchos::ScalarTraits<SC>::eps();
419 else
420 replaceWith = Teuchos::ScalarTraits<SC>::zero();
421
422 if(applyBCsTo22_) {
423 GetOStream(Runtime0) << "Maxwell1::compute(): nuking BC rows/cols of D0" << std::endl;
424 if (useKokkos_) {
425 Utilities_kokkos::ZeroDirichletCols(D0_Matrix_,BCcolsKokkos_,replaceWith);
426 } else {
427 Utilities::ZeroDirichletCols(D0_Matrix_,BCcols_,replaceWith);
428 }
429 }
430 else {
431 GetOStream(Runtime0) << "Maxwell1::compute(): nuking BC rows of D0" << std::endl;
432 if (useKokkos_) {
433 Utilities_kokkos::ZeroDirichletRows(D0_Matrix_,BCrowsKokkos_,replaceWith);
434 } else {
435 Utilities::ZeroDirichletRows(D0_Matrix_,BCrows_,replaceWith);
436 }
437 }
438
439 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
440 }
441
442
444 // What ML does is generate nodal prolongators with an auxillary hierarchy based on the
445 // user's (2,2) matrix. The actual nodal matrices for smoothing are generated by the
446 // Hiptmair smoother construction. We're not going to do that --- we'll
447 // do as we insert them into the final (1,1) hierarchy.
448
449 // Level 0
450 RCP<Matrix> Kn_Smoother_0;
451 if(have_generated_Kn) {
452 Kn_Smoother_0 = Kn_Matrix_;
453 }
454 else {
455 Kn_Smoother_0 = generate_kn();
456 }
457
459 // Copy the relevant (2,2) data to the (1,1) hierarchy
460 Hierarchy11_ = rcp(new Hierarchy("Maxwell1 (1,1)"));
461 RCP<Matrix> OldSmootherMatrix;
462 for(int i=0; i<Hierarchy22_->GetNumLevels(); i++) {
463 Hierarchy11_->AddNewLevel();
464 RCP<Level> NodeL = Hierarchy22_->GetLevel(i);
465 RCP<Level> EdgeL = Hierarchy11_->GetLevel(i);
466 RCP<Operator> NodeOp = NodeL->Get<RCP<Operator> >("A");
467 RCP<Matrix> NodeAggMatrix = rcp_dynamic_cast<Matrix>(NodeOp);
468 std::string labelstr = FormattingHelper::getColonLabel(EdgeL->getObjectLabel());
469
470 if(i==0) {
471 EdgeL->Set("A", SM_Matrix_);
472 EdgeL->Set("D0", D0_Matrix_);
473
474 EdgeL->Set("NodeAggMatrix",NodeAggMatrix);
475 EdgeL->Set("NodeMatrix",Kn_Smoother_0);
476 OldSmootherMatrix = Kn_Smoother_0;
477 }
478 else {
479 // Set the Nodal P
480 auto NodalP = NodeL->Get<RCP<Matrix> >("P");
481 EdgeL->Set("Pnodal",NodalP);
482
483 // If we repartition a processor away, a RCP<Operator> is stuck
484 // on the level instead of an RCP<Matrix>
485 if(!NodeAggMatrix.is_null()) {
486 EdgeL->Set("NodeAggMatrix",NodeAggMatrix);
487 if(!have_generated_Kn) {
488 // The user gave us a Kn, so we'll need to create the smoother matrix via RAP
489 RCP<Matrix> NewKn = Maxwell_Utils<SC,LO,GO,NO>::PtAPWrapper(OldSmootherMatrix,NodalP,parameterList_,labelstr);
490 EdgeL->Set("NodeMatrix",NewKn);
491 OldSmootherMatrix = NewKn;
492 }
493 else {
494 // The user didn't give us a Kn, so we aggregate and smooth with the same matrix
495 EdgeL->Set("NodeMatrix",NodeAggMatrix);
496 }
497 }
498 else {
499 // We've partitioned things away.
500 EdgeL->Set("NodeMatrix",NodeOp);
501 EdgeL->Set("NodeAggMatrix",NodeOp);
502 }
503 }
504
505 // Get the importer if we have one (for repartitioning)
506 // This will get used in ReitzingerPFactory
507 if(NodeL->IsAvailable("Importer")) {
508 auto importer = NodeL->Get<RCP<const Import> >("Importer");
509 EdgeL->Set("NodeImporter",importer);
510 }
511 }// end Hierarchy22 loop
512
513
515 // Generating the (1,1) Hierarchy
516 std::string fine_smoother = precList11_.get<std::string>("smoother: type");
517 {
518 SM_Matrix_->setObjectLabel("A(1,1)");
519 precList11_.set("coarse: max size",1);
520 precList11_.set("max levels",Hierarchy22_->GetNumLevels());
521
522 const bool refmaxwellCoarseSolve = (precList11_.get<std::string>("coarse: type",
523 MasterList::getDefault<std::string>("coarse: type")) == "RefMaxwell");
524 if (refmaxwellCoarseSolve) {
525 GetOStream(Runtime0) << "Maxwell1::compute(): Will set up RefMaxwell coarse solver" << std::endl;
526 precList11_.set("coarse: type", "none");
527 }
528
529 // Rip off non-serializable data before validation
530 Teuchos::ParameterList nonSerialList11, processedPrecList11;
531 MueLu::ExtractNonSerializableData(precList11_, processedPrecList11, nonSerialList11);
532 RCP<HierarchyManager<SC,LO,GO,NO> > mueLuFactory = rcp(new ParameterListInterpreter<SC,LO,GO,NO>(processedPrecList11,SM_Matrix_->getDomainMap()->getComm()));
533 Hierarchy11_->setlib(SM_Matrix_->getDomainMap()->lib());
534 Hierarchy11_->SetProcRankVerbose(SM_Matrix_->getDomainMap()->getComm()->getRank());
535 // Stick the non-serializible data on the hierarchy.
536 HierarchyUtils<SC,LO,GO,NO>::AddNonSerializableDataToHierarchy(*mueLuFactory,*Hierarchy11_, nonSerialList11);
537 mueLuFactory->SetupHierarchy(*Hierarchy11_);
538
539 if (refmaxwellCoarseSolve) {
540 GetOStream(Runtime0) << "Maxwell1::compute(): Setting up RefMaxwell coarse solver" << std::endl;
541 auto coarseLvl = Hierarchy11_->GetLevel(Hierarchy11_->GetNumLevels()-1);
542 auto coarseSolver = rcp(new MueLu::RefMaxwellSmoother<Scalar,LocalOrdinal,GlobalOrdinal,Node>("RefMaxwell",
543 precList11_.sublist("coarse: params")));
544 coarseSolver->Setup(*coarseLvl);
545 coarseLvl->Set("PreSmoother",
546 rcp_dynamic_cast<SmootherBase>(coarseSolver, true));
547 }
548
549 if(mode_ == MODE_REFMAXWELL) {
550 if(Hierarchy11_->GetNumLevels() > 1) {
551 RCP<Level> EdgeL = Hierarchy11_->GetLevel(1);
552 P11_ = EdgeL->Get<RCP<Matrix> >("P");
553 }
554 }
555 }
556
558 // Allocate temporary MultiVectors for solve (only needed for RefMaxwell style)
559 allocateMemory(1);
560
561 describe(GetOStream(Runtime0));
562
563
564#ifdef MUELU_MAXWELL1_DEBUG
565 for(int i=0; i<Hierarchy11_->GetNumLevels(); i++) {
566 RCP<Level> L = Hierarchy11_->GetLevel(i);
567 RCP<Matrix> EdgeMatrix = rcp_dynamic_cast<Matrix>(L->Get<RCP<Operator> >("A"));
568 RCP<Matrix> NodeMatrix = rcp_dynamic_cast<Matrix>(L->Get<RCP<Operator> >("NodeMatrix"));
569 RCP<Matrix> NodeAggMatrix = rcp_dynamic_cast<Matrix>(L->Get<RCP<Operator> >("NodeAggMatrix"));
570 RCP<Matrix> D0 =rcp_dynamic_cast<Matrix>( L->Get<RCP<Operator> >("D0"));
571
572 auto nrmE = EdgeMatrix->getFrobeniusNorm();
573 auto nrmN = NodeMatrix->getFrobeniusNorm();
574 auto nrmNa = NodeAggMatrix->getFrobeniusNorm();
575 auto nrmD0= D0->getFrobeniusNorm();
576
577 std::cout<<"DEBUG: Norms on Level "<<i<<" E/N/NA/D0 = "<<nrmE<<" / "<<nrmN <<" / "<<nrmNa<<" / "<< nrmD0 <<std::endl;
578 std::cout<<"DEBUG: NNZ on Level "<<i<<" E/N/NA/D0 = "<<
579 EdgeMatrix->getGlobalNumEntries()<<" / "<<
580 NodeMatrix->getGlobalNumEntries()<<" / "<<
581 NodeAggMatrix->getGlobalNumEntries()<<" / "<<
582 D0->getGlobalNumEntries()<<std::endl;
583 }
584#endif
585
586 }
587
588
589 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
590 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::generate_kn() const {
591 // NOTE: This does not nicely support reuse, but the relevant code can be copied from
592 // RefMaxwell when we decide we want to do this.
593
594 // NOTE: Boundary conditions OAZ are handled via the "rap: fix zero diagonals threshold"
595 RCP<Teuchos::TimeMonitor> tm = getTimer("MueLu Maxwell1: Build Kn");
596
597 Level fineLevel, coarseLevel;
598 fineLevel.SetFactoryManager(null);
599 coarseLevel.SetFactoryManager(null);
600 coarseLevel.SetPreviousLevel(rcpFromRef(fineLevel));
601 fineLevel.SetLevelID(0);
602 coarseLevel.SetLevelID(1);
603 fineLevel.Set("A",SM_Matrix_);
604 coarseLevel.Set("P",D0_Matrix_);
605 //coarseLevel.Set("Coordinates",Coords_);
606
607 coarseLevel.setlib(SM_Matrix_->getDomainMap()->lib());
608 fineLevel.setlib(SM_Matrix_->getDomainMap()->lib());
609 coarseLevel.setObjectLabel("Maxwell1 (2,2)");
610 fineLevel.setObjectLabel("Maxwell1 (2,2)");
611
612 RCP<RAPFactory> rapFact = rcp(new RAPFactory());
613 ParameterList rapList = *(rapFact->GetValidParameterList());
614 rapList.set("transpose: use implicit", true);
615 rapList.set("rap: triple product", parameterList_.get<bool>("rap: triple product", false));
616 rapFact->SetParameterList(rapList);
617 coarseLevel.Request("A", rapFact.get());
618 if (enable_reuse_) {
619 coarseLevel.Request("AP reuse data", rapFact.get());
620 coarseLevel.Request("RAP reuse data", rapFact.get());
621 }
622
623 RCP<Matrix> Kn_Matrix = coarseLevel.Get< RCP<Matrix> >("A", rapFact.get());
624 Kn_Matrix->setObjectLabel("A(2,2)");
625
626 return Kn_Matrix;
627 }
628
629
630
631 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
633 if(mode_ == MODE_REFMAXWELL) {
634 RCP<Teuchos::TimeMonitor> tmAlloc = getTimer("MueLu Maxwell1: Allocate MVs");
635
636 residualFine_ = MultiVectorFactory::Build(SM_Matrix_->getRangeMap(), numVectors);
637
638 if(!Hierarchy11_.is_null() && Hierarchy11_->GetNumLevels() > 1) {
639 RCP<Level> EdgeL = Hierarchy11_->GetLevel(1);
640 RCP<Matrix> A = EdgeL->Get<RCP<Matrix> >("A");
641 residual11c_ = MultiVectorFactory::Build(A->getRangeMap(), numVectors);
642 update11c_ = MultiVectorFactory::Build(A->getDomainMap(), numVectors);
643 }
644
645 if(!Hierarchy22_.is_null()) {
646 residual22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
647 update22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
648 }
649
650 }
651
652 }
653
654
655 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
656 void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dump(const Matrix& A, std::string name) const {
657 if (dump_matrices_) {
658 GetOStream(Runtime0) << "Dumping to " << name << std::endl;
659 Xpetra::IO<SC, LO, GO, NO>::Write(name, A);
660 }
661 }
662
663
664 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
665 void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dump(const MultiVector& X, std::string name) const {
666 if (dump_matrices_) {
667 GetOStream(Runtime0) << "Dumping to " << name << std::endl;
668 Xpetra::IO<SC, LO, GO, NO>::Write(name, X);
669 }
670 }
671
672
673 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
675 if (dump_matrices_) {
676 GetOStream(Runtime0) << "Dumping to " << name << std::endl;
677 Xpetra::IO<coordinateType, LO, GO, NO>::Write(name, X);
678 }
679 }
680
681
682 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
683 void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dump(const Teuchos::ArrayRCP<bool>& v, std::string name) const {
684 if (dump_matrices_) {
685 GetOStream(Runtime0) << "Dumping to " << name << std::endl;
686 std::ofstream out(name);
687 for (size_t i = 0; i < Teuchos::as<size_t>(v.size()); i++)
688 out << v[i] << "\n";
689 }
690 }
691
692 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
693 void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dump(const Kokkos::View<bool*, typename Node::device_type>& v, std::string name) const {
694 if (dump_matrices_) {
695 GetOStream(Runtime0) << "Dumping to " << name << std::endl;
696 std::ofstream out(name);
697 auto vH = Kokkos::create_mirror_view (v);
698 Kokkos::deep_copy(vH , v);
699 for (size_t i = 0; i < vH.size(); i++)
700 out << vH[i] << "\n";
701 }
702 }
703
704 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
705 Teuchos::RCP<Teuchos::TimeMonitor> Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getTimer(std::string name, RCP<const Teuchos::Comm<int> > comm) const {
706 if (IsPrint(Timings)) {
707 if (!syncTimers_)
708 return Teuchos::rcp(new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name)));
709 else {
710 if (comm.is_null())
711 return Teuchos::rcp(new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), SM_Matrix_->getRowMap()->getComm().ptr()));
712 else
713 return Teuchos::rcp(new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), comm.ptr()));
714 }
715 } else
716 return Teuchos::null;
717 }
718
719
720
721
722 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
723 void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::resetMatrix(RCP<Matrix> SM_Matrix_new, bool ComputePrec) {
724 bool reuse = !SM_Matrix_.is_null();
725 SM_Matrix_ = SM_Matrix_new;
726 dump(*SM_Matrix_, "SM.m");
727 if (ComputePrec)
728 compute(reuse);
729 }
730
731
732 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
734 // make sure that we have enough temporary memory
735 const SC one = Teuchos::ScalarTraits<SC>::one();
736 if (!allEdgesBoundary_ && X.getNumVectors() != residualFine_->getNumVectors())
737 allocateMemory(X.getNumVectors());
738
739 TEUCHOS_TEST_FOR_EXCEPTION(Hierarchy11_.is_null() || Hierarchy11_->GetNumLevels() == 0, Exceptions::RuntimeError, "(1,1) Hiearchy is null.");
740
741 // 1) Run fine pre-smoother using Hierarchy11
742 RCP<Level> Fine = Hierarchy11_->GetLevel(0);
743 if (Fine->IsAvailable("PreSmoother")) {
744 RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: PreSmoother");
745 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
746 preSmoo->Apply(X, RHS, true);
747 }
748
749 // 2) Compute residual
750 {
751 RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: residual calculation");
752 Utilities::Residual(*SM_Matrix_, X, RHS,*residualFine_);
753 }
754
755 // 3a) Restrict residual to (1,1) Hierarchy's level 1 and execute (1,1) hierarchy (use startLevel and InitialGuessIsZero)
756 if(!P11_.is_null()) {
757 RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: (1,1) correction");
758 P11_->apply(*residualFine_,*residual11c_,Teuchos::TRANS);
759 Hierarchy11_->Iterate(*residual11c_,*update11c_,true,1);
760 }
761
762 // 3b) Restrict residual to (2,2) Hierarchy's level 0 and execute (2,2) hierarchy (use InitialGuessIsZero)
763 if (!allNodesBoundary_) {
764 RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: (2,2) correction");
765 D0_Matrix_->apply(*residualFine_,*residual22_,Teuchos::TRANS);
766 Hierarchy22_->Iterate(*residual22_,*update22_,true,0);
767 }
768
769 // 4) Prolong both updates back into X-vector (Need to do both the P11 null and not null cases
770 {
771 RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: Orolongation");
772 if(!P11_.is_null())
773 P11_->apply(*update11c_,X,Teuchos::NO_TRANS,one,one);
774 if (!allNodesBoundary_)
775 D0_Matrix_->apply(*update22_,X,Teuchos::NO_TRANS,one,one);
776 }
777
778 // 5) Run fine post-smoother using Hierarchy11
779 if (Fine->IsAvailable("PostSmoother")) {
780 RCP<Teuchos::TimeMonitor> tmRes = getTimer("MueLu Maxwell1: PostSmoother");
781 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
782 postSmoo->Apply(X, RHS, false);
783 }
784
785
786 }
787
788 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
789 void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::applyInverseStandard(const MultiVector& RHS, MultiVector& X) const {
790 Hierarchy11_->Iterate(RHS,X,1,true);
791 }
792
793 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
794 void Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>::apply (const MultiVector& RHS, MultiVector& X,
795 Teuchos::ETransp /* mode */,
796 Scalar /* alpha */,
797 Scalar /* beta */) const {
798 RCP<Teuchos::TimeMonitor> tm = getTimer("MueLu Maxwell1: solve");
799 if(mode_ == MODE_GMHD_STANDARD)
800 HierarchyGmhd_->Iterate(RHS,X,1,true);
801 else if(mode_ == MODE_STANDARD || mode_ == MODE_EDGE_ONLY)
802 applyInverseStandard(RHS,X);
803 else if(mode_ == MODE_REFMAXWELL)
804 applyInverseRefMaxwellAdditive(RHS,X);
805 else
806 TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Must use mode 'standard', 'refmaxwell' or 'edge only' when not doing GMHD.");
807 }
808
809
810
811 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
813 return false;
814 }
815
816
817 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
819 initialize(const Teuchos::RCP<Matrix> & D0_Matrix,
820 const Teuchos::RCP<Matrix> & Kn_Matrix,
821 const Teuchos::RCP<MultiVector> & Nullspace,
822 const Teuchos::RCP<RealValuedMultiVector> & Coords,
823 Teuchos::ParameterList& List)
824 {
825 // some pre-conditions
826 TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
827
828#ifdef HAVE_MUELU_DEBUG
829 if(!Kn_Matrix.is_null()) {
830 TEUCHOS_ASSERT(Kn_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getDomainMap()));
831 TEUCHOS_ASSERT(Kn_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getDomainMap()));
832 }
833
834
835 TEUCHOS_ASSERT(D0_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getRowMap()));
836#endif
837
838 Hierarchy11_ = Teuchos::null;
839 Hierarchy22_ = Teuchos::null;
840 HierarchyGmhd_ = Teuchos::null;
841 if (mode_ != MODE_GMHD_STANDARD) mode_ = MODE_STANDARD;
842
843 // Default settings
844 useKokkos_=false;
845 allEdgesBoundary_=false;
846 allNodesBoundary_=false;
847 dump_matrices_ = false;
848 enable_reuse_=false;
849 syncTimers_=false;
850 applyBCsTo22_ = false;
851
852
853 // set parameters
854 setParameters(List);
855
856 if (D0_Matrix->getRowMap()->lib() == Xpetra::UseTpetra) {
857 // We will remove boundary conditions from D0, and potentially change maps, so we copy the input.
858 // Fortunately, D0 is quite sparse.
859 // We cannot use the Tpetra copy constructor, since it does not copy the graph.
860
861 RCP<Matrix> D0copy = MatrixFactory::Build(D0_Matrix->getRowMap(), D0_Matrix->getColMap(), 0);
862 RCP<CrsMatrix> D0copyCrs = rcp_dynamic_cast<CrsMatrixWrap>(D0copy,true)->getCrsMatrix();
863 ArrayRCP<const size_t> D0rowptr_RCP;
864 ArrayRCP<const LO> D0colind_RCP;
865 ArrayRCP<const SC> D0vals_RCP;
866 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,true)->getCrsMatrix()->getAllValues(D0rowptr_RCP,
867 D0colind_RCP,
868 D0vals_RCP);
869
870 ArrayRCP<size_t> D0copyrowptr_RCP;
871 ArrayRCP<LO> D0copycolind_RCP;
872 ArrayRCP<SC> D0copyvals_RCP;
873 D0copyCrs->allocateAllValues(D0vals_RCP.size(),D0copyrowptr_RCP,D0copycolind_RCP,D0copyvals_RCP);
874 D0copyrowptr_RCP.deepCopy(D0rowptr_RCP());
875 D0copycolind_RCP.deepCopy(D0colind_RCP());
876 D0copyvals_RCP.deepCopy(D0vals_RCP());
877 D0copyCrs->setAllValues(D0copyrowptr_RCP,
878 D0copycolind_RCP,
879 D0copyvals_RCP);
880 D0copyCrs->expertStaticFillComplete(D0_Matrix->getDomainMap(), D0_Matrix->getRangeMap(),
881 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,true)->getCrsMatrix()->getCrsGraph()->getImporter(),
882 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,true)->getCrsMatrix()->getCrsGraph()->getExporter());
883 D0_Matrix_ = D0copy;
884 } else
885 D0_Matrix_ = MatrixFactory::BuildCopy(D0_Matrix);
886
887
888 Kn_Matrix_ = Kn_Matrix;
889 Coords_ = Coords;
890 Nullspace_ = Nullspace;
891
892 if(!Kn_Matrix_.is_null()) dump(*Kn_Matrix_, "Kn.m");
893 if (!Nullspace_.is_null()) dump(*Nullspace_, "nullspace.m");
894 if (!Coords_.is_null()) dumpCoords(*Coords_, "coords.m");
895
896 }
897
898
899
900 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
902 describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel /* verbLevel */) const {
903
904 std::ostringstream oss;
905
906 RCP<const Teuchos::Comm<int> > comm = SM_Matrix_->getDomainMap()->getComm();
907
908#ifdef HAVE_MPI
909 int root;
910 if (!Kn_Matrix_.is_null())
911 root = comm->getRank();
912 else
913 root = -1;
914
915 int actualRoot;
916 reduceAll(*comm, Teuchos::REDUCE_MAX, root, Teuchos::ptr(&actualRoot));
917 root = actualRoot;
918#endif
919
920
921 oss << "\n--------------------------------------------------------------------------------\n" <<
922 "--- Maxwell1 Summary ---\n"
923 "--------------------------------------------------------------------------------" << std::endl;
924 oss << std::endl;
925
926 GlobalOrdinal numRows;
927 GlobalOrdinal nnz;
928
929 SM_Matrix_->getRowMap()->getComm()->barrier();
930
931 numRows = SM_Matrix_->getGlobalNumRows();
932 nnz = SM_Matrix_->getGlobalNumEntries();
933
934 Xpetra::global_size_t tt = numRows;
935 int rowspacer = 3; while (tt != 0) { tt /= 10; rowspacer++; }
936 tt = nnz;
937 int nnzspacer = 2; while (tt != 0) { tt /= 10; nnzspacer++; }
938
939 oss << "block " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz " << std::setw(9) << " nnz/row" << std::endl;
940 oss << "(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
941
942 if (!Kn_Matrix_.is_null()) {
943 // ToDo: make sure that this is printed correctly
944 numRows = Kn_Matrix_->getGlobalNumRows();
945 nnz = Kn_Matrix_->getGlobalNumEntries();
946
947 oss << "(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
948 }
949
950 oss << std::endl;
951
952
953 std::string outstr = oss.str();
954
955#ifdef HAVE_MPI
956 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
957 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
958
959 int strLength = outstr.size();
960 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
961 if (comm->getRank() != root)
962 outstr.resize(strLength);
963 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
964#endif
965
966 out << outstr;
967
968 if (!Hierarchy11_.is_null())
969 Hierarchy11_->describe(out, GetVerbLevel());
970
971 if (!Hierarchy22_.is_null())
972 Hierarchy22_->describe(out, GetVerbLevel());
973
974 if (!HierarchyGmhd_.is_null())
975 HierarchyGmhd_->describe(out, GetVerbLevel());
976
977 if (IsPrint(Statistics2)) {
978 // Print the grid of processors
979 std::ostringstream oss2;
980
981 oss2 << "Sub-solver distribution over ranks" << std::endl;
982 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;
983
984 int numProcs = comm->getSize();
985#ifdef HAVE_MPI
986 RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
987
988 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
989#endif
990
991 char status = 0;
992 if (!Kn_Matrix_.is_null())
993 status += 1;
994 std::vector<char> states(numProcs, 0);
995#ifdef HAVE_MPI
996 MPI_Gather(&status, 1, MPI_CHAR, &states[0], 1, MPI_CHAR, 0, *rawMpiComm);
997#else
998 states.push_back(status);
999#endif
1000
1001 int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
1002 for (int proc = 0; proc < numProcs; proc += rowWidth) {
1003 for (int j = 0; j < rowWidth; j++)
1004 if (proc + j < numProcs)
1005 if (states[proc+j] == 0)
1006 oss2 << ".";
1007 else if (states[proc+j] == 1)
1008 oss2 << "1";
1009 else if (states[proc+j] == 2)
1010 oss2 << "2";
1011 else
1012 oss2 << "B";
1013 else
1014 oss2 << " ";
1015
1016 oss2 << " " << proc << ":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
1017 }
1018 oss2 << std::endl;
1019 GetOStream(Statistics2) << oss2.str();
1020 }
1021
1022
1023 }
1024
1025
1026
1027} // namespace
1028
1029#define MUELU_MAXWELL1_SHORT
1030#endif //ifdef MUELU_MAXWELL1_DEF_HPP
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
Exception throws to report errors in the internal logical of the program.
static void AddNonSerializableDataToHierarchy(HierarchyManager &HM, Hierarchy &H, const ParameterList &nonSerialList)
Add non-serializable data to Hierarchy.
static void CopyBetweenHierarchies(Hierarchy &fromHierarchy, Hierarchy &toHierarchy, const std::string fromLabel, const std::string toLabel, const std::string dataType)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void setlib(Xpetra::UnderlyingLib lib2)
void SetLevelID(int levelID)
Set level number.
Definition: MueLu_Level.cpp:78
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
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
static std::string translate(Teuchos::ParameterList &paramList, const std::string &defaultVals="")
: Translate ML parameters to MueLu parameter XML string
Xpetra::MultiVector< coordinateType, LO, GO, NO > RealValuedMultiVector
Teuchos::RCP< Matrix > generate_kn() const
Generates the Kn matrix.
Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
void compute(bool reuse=false)
Setup the preconditioner.
void GMHDSetupHierarchy(Teuchos::ParameterList &List) const
Sets up hiearchy for GMHD matrices that include generalized Ohms law equations.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
void applyInverseRefMaxwellAdditive(const MultiVector &RHS, MultiVector &X) const
apply RefMaxwell additive 2x2 style cycle
void applyInverseStandard(const MultiVector &RHS, MultiVector &X) const
apply standard Maxwell1 cycle
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)
void dumpCoords(const RealValuedMultiVector &X, std::string name) const
dump out real-valued multivector
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
void dump(const Matrix &A, std::string name) const
dump out matrix
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
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
void setParameters(Teuchos::ParameterList &list)
Set parameters.
void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
void allocateMemory(int numVectors) const
allocate multivectors for solve
static void detectBoundaryConditionsSM(RCP< Matrix > &SM_Matrix, RCP< Matrix > &D0_Matrix, magnitudeType rowSumTol, bool useKokkos_, Kokkos::View< bool *, typename Node::device_type > &BCrowsKokkos, Kokkos::View< bool *, typename Node::device_type > &BCcolsKokkos, Kokkos::View< bool *, typename Node::device_type > &BCdomainKokkos, int &BCedges, int &BCnodes, Teuchos::ArrayRCP< bool > &BCrows, Teuchos::ArrayRCP< bool > &BCcols, Teuchos::ArrayRCP< bool > &BCdomain, bool &allEdgesBoundary, bool &allNodesBoundary)
Detect Dirichlet boundary conditions.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > PtAPWrapper(RCP< Matrix > &A, RCP< Matrix > &P, Teuchos::ParameterList &params, std::string &label)
Factory for building coarse matrices.
Class that encapsulates Operator smoothers.
static void ZeroDirichletRows(RCP< Matrix > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static void ZeroDirichletCols(RCP< Matrix > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletCols, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static void ZeroDirichletCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Xpetra::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)
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
Namespace for MueLu classes and methods.
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.
@ Warnings0
Important warning messages (one line)
@ Statistics2
Print even more statistics.
@ Runtime0
One-liner description of what is happening.
@ Warnings1
Additional warnings.
@ Timings
Print all timing information.
MsgType toVerbLevel(const std::string &verbLevelStr)
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 std::string getColonLabel(const std::string &label)
Helper function for object label.