MueLu  Version of the Day
MueLu_Maxwell_Utils_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_MAXWELL_UTILS_DEF_HPP
47 #define MUELU_MAXWELL_UTILS_DEF_HPP
48 
49 #include <sstream>
50 
51 #include "MueLu_ConfigDefs.hpp"
52 
53 #include "Xpetra_Map.hpp"
54 #include "Xpetra_CrsMatrixUtils.hpp"
55 #include "Xpetra_MatrixUtils.hpp"
56 
58 
59 #include "MueLu_Level.hpp"
60 #include "MueLu_ThresholdAFilterFactory.hpp"
61 #include "MueLu_Utilities.hpp"
62 
63 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
64 #include "MueLu_Utilities_kokkos.hpp"
65 #endif
66 
67 // Stratimikos
68 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
69 // Thyra includes
70 #include <Thyra_VectorBase.hpp>
71 #include <Thyra_SolveSupportTypes.hpp>
72 // Stratimikos includes
73 #include <Stratimikos_DefaultLinearSolverBuilder.hpp>
75 #include "Teuchos_AbstractFactoryStd.hpp"
76 // Ifpack2 includes
77 #ifdef HAVE_MUELU_IFPACK2
78 #include <Thyra_Ifpack2PreconditionerFactory.hpp>
79 #endif
80 #endif
81 
82 namespace MueLu {
83 
84 
85  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
87  RCP<Matrix> & D0_Matrix_,
88  magnitudeType rowSumTol,
90  bool useKokkos_,
94 #endif
95  int & BCedges_,
96  int & BCnodes_,
97  Teuchos::ArrayRCP<bool> & BCrows_,
98  Teuchos::ArrayRCP<bool> & BCcols_,
99  Teuchos::ArrayRCP<bool> & BCdomain_,
100  bool & allEdgesBoundary_,
101  bool & allNodesBoundary_) {
102  // clean rows associated with boundary conditions
103  // Find rows with only 1 or 2 nonzero entries, record them in BCrows_.
104  // BCrows_[i] is true, iff i is a boundary row
105  // BCcols_[i] is true, iff i is a boundary column
106  int BCedgesLocal = 0;
107  int BCnodesLocal = 0;
108 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
109  if (useKokkos_) {
110  BCrowsKokkos_ = Utilities_kokkos::DetectDirichletRows(*SM_Matrix_,Teuchos::ScalarTraits<magnitudeType>::eps(),/*count_twos_as_dirichlet=*/true);
111 
112  if (rowSumTol > 0.)
113  Utilities_kokkos::ApplyRowSumCriterion(*SM_Matrix_, rowSumTol, BCrowsKokkos_);
114 
115  BCcolsKokkos_ = Kokkos::View<bool*,typename Node::device_type>(Kokkos::ViewAllocateWithoutInitializing("dirichletCols"), D0_Matrix_->getColMap()->getNodeNumElements());
116  BCdomainKokkos_ = Kokkos::View<bool*,typename Node::device_type>(Kokkos::ViewAllocateWithoutInitializing("dirichletDomains"), D0_Matrix_->getDomainMap()->getNodeNumElements());
117  Utilities_kokkos::DetectDirichletColsAndDomains(*D0_Matrix_,BCrowsKokkos_,BCcolsKokkos_,BCdomainKokkos_);
118 
119  auto BCrowsKokkos=BCrowsKokkos_;
120  Kokkos::parallel_reduce(BCrowsKokkos_.size(), KOKKOS_LAMBDA (int i, int & sum) {
121  if (BCrowsKokkos(i))
122  ++sum;
123  }, BCedgesLocal );
124 
125  auto BCdomainKokkos = BCdomainKokkos_;
126  Kokkos::parallel_reduce(BCdomainKokkos_.size(), KOKKOS_LAMBDA (int i, int & sum) {
127  if (BCdomainKokkos(i))
128  ++sum;
129  }, BCnodesLocal);
130  } else
131 #endif // HAVE_MUELU_KOKKOS_REFACTOR
132  {
133  BCrows_ = Teuchos::arcp_const_cast<bool>(Utilities::DetectDirichletRows(*SM_Matrix_,Teuchos::ScalarTraits<magnitudeType>::eps(),/*count_twos_as_dirichlet=*/true));
134 
135  if (rowSumTol > 0.)
136  Utilities::ApplyRowSumCriterion(*SM_Matrix_, rowSumTol, BCrows_);
137 
138  BCcols_.resize(D0_Matrix_->getColMap()->getNodeNumElements());
139  BCdomain_.resize(D0_Matrix_->getDomainMap()->getNodeNumElements());
140  Utilities::DetectDirichletColsAndDomains(*D0_Matrix_,BCrows_,BCcols_,BCdomain_);
141 
142  for (auto it = BCrows_.begin(); it != BCrows_.end(); ++it)
143  if (*it)
144  BCedgesLocal += 1;
145  for (auto it = BCdomain_.begin(); it != BCdomain_.end(); ++it)
146  if (*it)
147  BCnodesLocal += 1;
148  }
149 
150  MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCedgesLocal, BCedges_);
151  MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCnodesLocal, BCnodes_);
152 
153 
154  allEdgesBoundary_ = Teuchos::as<Xpetra::global_size_t>(BCedges_) >= D0_Matrix_->getRangeMap()->getGlobalNumElements();
155  allNodesBoundary_ = Teuchos::as<Xpetra::global_size_t>(BCnodes_) >= D0_Matrix_->getDomainMap()->getGlobalNumElements();
156  }
157 
158 
159  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
161  RCP<Matrix> & D0_Matrix_,
162  RCP<Matrix> & SM_Matrix_,
163  RCP<Matrix> & M1_Matrix_,
164  RCP<Matrix> & Ms_Matrix_) {
165 
166  bool defaultFilter = false;
167 
168  // Remove zero entries from D0 if necessary.
169  // In the construction of the prolongator we use the graph of the
170  // matrix, so zero entries mess it up.
171  if (parameterList_.get<bool>("refmaxwell: filter D0", true) && D0_Matrix_->getNodeMaxNumRowEntries()>2) {
172  Level fineLevel;
173  fineLevel.SetFactoryManager(null);
174  fineLevel.SetLevelID(0);
175  fineLevel.Set("A",D0_Matrix_);
176  fineLevel.setlib(D0_Matrix_->getDomainMap()->lib());
177  // We expect D0 to have entries +-1, so any threshold value will do.
178  RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",1.0e-8,/*keepDiagonal=*/false,/*expectedNNZperRow=*/2));
179  fineLevel.Request("A",ThreshFact.get());
180  ThreshFact->Build(fineLevel);
181  D0_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
182 
183  // If D0 has too many zeros, maybe SM and M1 do as well.
184  defaultFilter = true;
185  }
186 
187  if (!M1_Matrix_.is_null() && parameterList_.get<bool>("refmaxwell: filter M1", defaultFilter)) {
188  RCP<Vector> diag = VectorFactory::Build(M1_Matrix_->getRowMap());
189  // find a reasonable absolute value threshold
190  M1_Matrix_->getLocalDiagCopy(*diag);
191  magnitudeType threshold = 1.0e-8 * diag->normInf();
192 
193  Level fineLevel;
194  fineLevel.SetFactoryManager(null);
195  fineLevel.SetLevelID(0);
196  fineLevel.Set("A",M1_Matrix_);
197  fineLevel.setlib(M1_Matrix_->getDomainMap()->lib());
198  RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",threshold,/*keepDiagonal=*/true));
199  fineLevel.Request("A",ThreshFact.get());
200  ThreshFact->Build(fineLevel);
201  M1_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
202  }
203 
204  if (!Ms_Matrix_.is_null() && parameterList_.get<bool>("refmaxwell: filter Ms", defaultFilter)) {
205  RCP<Vector> diag = VectorFactory::Build(Ms_Matrix_->getRowMap());
206  // find a reasonable absolute value threshold
207  Ms_Matrix_->getLocalDiagCopy(*diag);
208  magnitudeType threshold = 1.0e-8 * diag->normInf();
209 
210  Level fineLevel;
211  fineLevel.SetFactoryManager(null);
212  fineLevel.SetLevelID(0);
213  fineLevel.Set("A",Ms_Matrix_);
214  fineLevel.setlib(Ms_Matrix_->getDomainMap()->lib());
215  RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",threshold,/*keepDiagonal=*/true));
216  fineLevel.Request("A",ThreshFact.get());
217  ThreshFact->Build(fineLevel);
218  Ms_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
219  }
220 
221  if (!SM_Matrix_.is_null() && parameterList_.get<bool>("refmaxwell: filter SM", defaultFilter)) {
222  RCP<Vector> diag = VectorFactory::Build(SM_Matrix_->getRowMap());
223  // find a reasonable absolute value threshold
224  SM_Matrix_->getLocalDiagCopy(*diag);
225  magnitudeType threshold = 1.0e-8 * diag->normInf();
226 
227  Level fineLevel;
228  fineLevel.SetFactoryManager(null);
229  fineLevel.SetLevelID(0);
230  fineLevel.Set("A",SM_Matrix_);
231  fineLevel.setlib(SM_Matrix_->getDomainMap()->lib());
232  RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",threshold,/*keepDiagonal=*/true));
233  fineLevel.Request("A",ThreshFact.get());
234  ThreshFact->Build(fineLevel);
235  SM_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
236  }
237 
238  }
239 
240 
241 
242  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
244  setMatvecParams(Matrix& A, RCP<ParameterList> matvecParams) {
245  RCP<const Import> xpImporter = A.getCrsGraph()->getImporter();
246  if (!xpImporter.is_null())
247  xpImporter->setDistributorParameters(matvecParams);
248  RCP<const Export> xpExporter = A.getCrsGraph()->getExporter();
249  if (!xpExporter.is_null())
250  xpExporter->setDistributorParameters(matvecParams);
251  }
252 
253 
254 
255 
256 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
257 
258  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
259  struct StratimikosWrapperImpl {
260  static RCP<Thyra::PreconditionerBase<Scalar> > setupStratimikosPreconditioner(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A,
261  RCP<ParameterList> params) {
262  throw std::runtime_error("setupStratimikosPreconditioner: Requires Scalar=double");
263  }
264  };
265 
266  template<class LocalOrdinal, class GlobalOrdinal, class Node>
267  struct StratimikosWrapperImpl<double,LocalOrdinal,GlobalOrdinal,Node> {
268  static RCP<Thyra::PreconditionerBase<double> > setupStratimikosPreconditioner(RCP<Xpetra::Matrix<double,LocalOrdinal,GlobalOrdinal,Node> > A,
269  RCP<ParameterList> params) {
270  typedef double Scalar;
271 
272  // Build Thyra linear algebra objects
273  RCP<const Thyra::LinearOpBase<Scalar> > thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(A)->getCrsMatrix());
274 
275  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
276  typedef Thyra::PreconditionerFactoryBase<Scalar> Base;
277  typedef Thyra::MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> ImplMueLu;
278  linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, ImplMueLu>(), "MueLu");
279 #ifdef HAVE_MUELU_IFPACK2
280  // Register Ifpack2 as a Stratimikos preconditioner strategy.
281  typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Impl;
282  linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(), "Ifpack2");
283 #endif
284 
285  linearSolverBuilder.setParameterList(params);
286 
287  // Build a new "solver factory" according to the previously specified parameter list.
288  // RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
289 
290  auto precFactory = Thyra::createPreconditioningStrategy(linearSolverBuilder);
291  auto prec = precFactory->createPrec();
292 
293  precFactory->initializePrec(Thyra::defaultLinearOpSource(thyraA), prec.get(), Thyra::SUPPORT_SOLVE_UNSPECIFIED);
294 
295  return prec;
296  }
297  };
298 
299 
300  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
301  RCP<Thyra::PreconditionerBase<Scalar> >
302  Maxwell_Utils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setupStratimikosPreconditioner(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > A,
303  RCP<ParameterList> params) {
304  return StratimikosWrapperImpl<SC,LO,GO,NO>::setupStratimikosPreconditioner(A,params);
305  }
306 
307 
308 
309 
310 #endif
311 
312 } // namespace
313 
314 #define MUELU_MAXWELL_UTILS_SHORT
315 #endif //ifdef MUELU_MAXWELL_UTILS_DEF_HPP
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
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 ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
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.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
Namespace for MueLu classes and methods.
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
Definition: MueLu_Level.cpp:92
MueLu::DefaultNode Node
void setlib(Xpetra::UnderlyingLib lib2)
void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Kokkos::View< bool *, typename Node::device_type > &dirichletRows, Kokkos::View< bool *, typename Node::device_type > dirichletCols, Kokkos::View< bool *, typename Node::device_type > dirichletDomain)
Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
static void setMatvecParams(Matrix &A, RCP< ParameterList > matvecParams)
Sets matvec params on a matrix.
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.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomains)
void SetLevelID(int levelID)
Set level number.
Definition: MueLu_Level.cpp:78
void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, typename Node::device_type > &dirichletRows)
Factory for building a thresholded operator.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.