MueLu  Version of the Day
MueLu_RefMaxwell_decl.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_REFMAXWELL_DECL_HPP
47 #define MUELU_REFMAXWELL_DECL_HPP
48 
49 #include "MueLu_ConfigDefs.hpp"
50 #include "MueLu_BaseClass.hpp"
51 
57 #include "MueLu_SaPFactory_fwd.hpp"
60 #include "MueLu_Utilities_fwd.hpp"
61 
62 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
70 #endif
71 
78 
80 #include "MueLu_TrilinosSmoother.hpp"
81 #include "MueLu_Hierarchy.hpp"
82 
83 #include "Xpetra_Map_fwd.hpp"
84 #include "Xpetra_Matrix_fwd.hpp"
85 #include "Xpetra_MatrixFactory_fwd.hpp"
86 #include "Xpetra_MultiVectorFactory_fwd.hpp"
87 #include "Xpetra_VectorFactory_fwd.hpp"
88 #include "Xpetra_CrsMatrixWrap_fwd.hpp"
89 
90 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA))
91 #define MUELU_REFMAXWELL_CAN_USE_HIPTMAIR
93 #include "Ifpack2_Hiptmair.hpp"
94 #endif
95 
96 // Stratimikos
97 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
98 #include <Thyra_LinearOpWithSolveBase.hpp>
99 #endif
100 
101 namespace MueLu {
102 
129  template <class Scalar,
130  class LocalOrdinal,
131  class GlobalOrdinal,
132  class Node>
133  class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
134 
135 #undef MUELU_REFMAXWELL_SHORT
136 #include "MueLu_UseShortNames.hpp"
137 
138  public:
139 
140  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
141  typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType coordinateType;
142  typedef typename Xpetra::MultiVector<coordinateType,LO,GO,NO> RealValuedMultiVector;
143 
146  HierarchyH_(Teuchos::null),
147  Hierarchy22_(Teuchos::null),
148  disable_addon_(MasterList::getDefault<bool>("refmaxwell: disable addon")),
149  mode_(MasterList::getDefault<std::string>("refmaxwell: mode"))
150  {
151  }
152 
154  RefMaxwell(Teuchos::RCP<Hierarchy> HH, Teuchos::RCP<Hierarchy> H22) :
155  HierarchyH_(HH),
156  Hierarchy22_(H22),
157  disable_addon_(MasterList::getDefault<bool>("refmaxwell: disable addon")),
158  mode_(MasterList::getDefault<std::string>("refmaxwell: mode"))
159  {
160  }
161 
174  RefMaxwell(const Teuchos::RCP<Matrix> & SM_Matrix,
175  const Teuchos::RCP<Matrix> & D0_Matrix,
176  const Teuchos::RCP<Matrix> & Ms_Matrix,
177  const Teuchos::RCP<Matrix> & M0inv_Matrix,
178  const Teuchos::RCP<Matrix> & M1_Matrix,
179  const Teuchos::RCP<MultiVector> & Nullspace,
180  const Teuchos::RCP<RealValuedMultiVector> & Coords,
181  Teuchos::ParameterList& List,
182  bool ComputePrec = true)
183  {
184  initialize(D0_Matrix,Ms_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
185  resetMatrix(SM_Matrix,ComputePrec);
186  }
187 
199  RefMaxwell(const Teuchos::RCP<Matrix> & SM_Matrix,
200  const Teuchos::RCP<Matrix> & D0_Matrix,
201  const Teuchos::RCP<Matrix> & M0inv_Matrix,
202  const Teuchos::RCP<Matrix> & M1_Matrix,
203  const Teuchos::RCP<MultiVector> & Nullspace,
204  const Teuchos::RCP<RealValuedMultiVector> & Coords,
205  Teuchos::ParameterList& List,
206  bool ComputePrec = true)
207  {
208  initialize(D0_Matrix,M1_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
209  resetMatrix(SM_Matrix,ComputePrec);
210  }
211 
221  RefMaxwell(const Teuchos::RCP<Matrix> & D0_Matrix,
222  const Teuchos::RCP<Matrix> & M0inv_Matrix,
223  const Teuchos::RCP<Matrix> & M1_Matrix,
224  const Teuchos::RCP<MultiVector> & Nullspace,
225  const Teuchos::RCP<RealValuedMultiVector> & Coords,
226  Teuchos::ParameterList& List) : SM_Matrix_(Teuchos::null)
227  {
228  initialize(D0_Matrix,M1_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
229  }
230 
241  RefMaxwell(const Teuchos::RCP<Matrix> & SM_Matrix,
242  const Teuchos::RCP<Matrix> & D0_Matrix,
243  const Teuchos::RCP<Matrix> & M1_Matrix,
244  const Teuchos::RCP<MultiVector> & Nullspace,
245  const Teuchos::RCP<RealValuedMultiVector> & Coords,
246  Teuchos::ParameterList& List,
247  bool ComputePrec)
248  {
249  initialize(D0_Matrix,M1_Matrix,Teuchos::null,M1_Matrix,Nullspace,Coords,List);
250  resetMatrix(SM_Matrix,ComputePrec);
251  }
252 
261  RefMaxwell(const Teuchos::RCP<Matrix> & D0_Matrix,
262  const Teuchos::RCP<Matrix> & M1_Matrix,
263  const Teuchos::RCP<MultiVector> & Nullspace,
264  const Teuchos::RCP<RealValuedMultiVector> & Coords,
265  Teuchos::ParameterList& List) : SM_Matrix_(Teuchos::null)
266  {
267  initialize(D0_Matrix,M1_Matrix,Teuchos::null,M1_Matrix,Nullspace,Coords,List);
268  }
269 
276  RefMaxwell(const Teuchos::RCP<Matrix> & SM_Matrix,
277  Teuchos::ParameterList& List,
278  bool ComputePrec = true)
279  {
280 
281  RCP<MultiVector> Nullspace = List.get<RCP<MultiVector> >("Nullspace", Teuchos::null);
282  RCP<RealValuedMultiVector> Coords = List.get<RCP<RealValuedMultiVector> >("Coordinates", Teuchos::null);
283  RCP<Matrix> D0_Matrix = List.get<RCP<Matrix> >("D0");
284  RCP<Matrix> Ms_Matrix;
285  if (List.isType<RCP<Matrix> >("Ms"))
286  Ms_Matrix = List.get<RCP<Matrix> >("Ms");
287  else
288  Ms_Matrix = List.get<RCP<Matrix> >("M1");
289  RCP<Matrix> M1_Matrix = List.get<RCP<Matrix> >("M1");
290  RCP<Matrix> M0inv_Matrix = List.get<RCP<Matrix> >("M0inv", Teuchos::null);
291 
292  initialize(D0_Matrix,Ms_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
293 
294  if (SM_Matrix != Teuchos::null)
295  resetMatrix(SM_Matrix,ComputePrec);
296  }
297 
299  virtual ~RefMaxwell() {}
300 
302  Teuchos::RCP<const Map> getDomainMap() const;
303 
305  Teuchos::RCP<const Map> getRangeMap() const;
306 
308  const Teuchos::RCP<Matrix> & getJacobian() const {
309  return SM_Matrix_;
310  }
311 
313  void setParameters(Teuchos::ParameterList& list);
314 
316  void compute(bool reuse=false);
317 
319  void buildProlongator();
320 
322  void formCoarseMatrix();
323 
325  void resetMatrix(Teuchos::RCP<Matrix> SM_Matrix_new, bool ComputePrec=true);
326 
330  void apply (const MultiVector& X, MultiVector& Y,
331  Teuchos::ETransp mode = Teuchos::NO_TRANS,
332  Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
333  Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const;
334 
336  bool hasTransposeApply() const;
337 
338  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_HIGH) const;
339 
341  void residual(const MultiVector & X,
342  const MultiVector & B,
343  MultiVector & R) const {
344  using STS = Teuchos::ScalarTraits<Scalar>;
345  R.update(STS::one(),B,STS::zero());
346  this->apply (X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
347  }
348 
349 
350  private:
351 
362  void initialize(const Teuchos::RCP<Matrix> & D0_Matrix,
363  const Teuchos::RCP<Matrix> & Ms_Matrix,
364  const Teuchos::RCP<Matrix> & M0inv_Matrix,
365  const Teuchos::RCP<Matrix> & M1_Matrix,
366  const Teuchos::RCP<MultiVector> & Nullspace,
367  const Teuchos::RCP<RealValuedMultiVector> & Coords,
368  Teuchos::ParameterList& List);
369 
371  void setFineLevelSmoother();
372 
374  void applyInverseAdditive(const MultiVector& RHS, MultiVector& X) const;
375 
377  void solveH(const MultiVector& RHS, MultiVector& X) const;
378 
380  void solve22(const MultiVector& RHS, MultiVector& X) const;
381 
383  void allocateMemory(int numVectors) const;
384 
386  void dump(const Matrix& A, std::string name) const;
387 
389  void dump(const MultiVector& X, std::string name) const;
390 
392  void dumpCoords(const RealValuedMultiVector& X, std::string name) const;
393 
395  void dump(const Teuchos::ArrayRCP<bool>& v, std::string name) const;
396 
397 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
398  void dump(const Kokkos::View<bool*, typename Node::device_type>& v, std::string name) const;
400 #endif
401 
403  Teuchos::RCP<Teuchos::TimeMonitor> getTimer(std::string name, RCP<const Teuchos::Comm<int> > comm=Teuchos::null) const;
404 
405 
407  Teuchos::RCP<Hierarchy> HierarchyH_, Hierarchy22_;
408  Teuchos::RCP<SmootherBase> PreSmoother_, PostSmoother_;
409  Teuchos::RCP<SmootherPrototype> PreSmootherData_, PostSmootherData_;
410 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR)
411  Teuchos::RCP<Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> > hiptmairPreSmoother_, hiptmairPostSmoother_;
412 #endif
414 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
415  RCP<Thyra::PreconditionerBase<Scalar> > thyraPrecH_, thyraPrec22_;
416 #endif
419  Teuchos::RCP<Matrix> A_nodal_Matrix_, P11_, R11_, AH_, A22_, Addon_Matrix_;
420  Teuchos::RCP<const Map> D0origDomainMap_;
421  Teuchos::RCP<const Import> D0origImporter_;
423 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
424  Kokkos::View<bool*, typename Node::device_type> BCrowsKokkos_, BCcolsKokkos_, BCdomainKokkos_;
425 #endif
427  Teuchos::ArrayRCP<bool> BCrows_, BCcols_, BCdomain_;
429  Teuchos::RCP<MultiVector> Nullspace_;
431  Teuchos::RCP<RealValuedMultiVector> Coords_, CoordsH_;
433  Teuchos::RCP<MultiVector> NullspaceH_;
435  Teuchos::RCP<const Import> ImporterH_, Importer22_;
440  Teuchos::RCP<Teuchos::ParameterList> AH_AP_reuse_data_, AH_RAP_reuse_data_;
441  Teuchos::RCP<Teuchos::ParameterList> A22_AP_reuse_data_, A22_RAP_reuse_data_;
446  std::string mode_;
449  };
450 
451 
452 } // namespace
453 
454 #define MUELU_REFMAXWELL_SHORT
455 #endif // MUELU_REFMAXWELL_DECL_HPP
Teuchos::RCP< SmootherBase > PostSmoother_
Teuchos::RCP< const Import > Importer22_
Teuchos::RCP< Matrix > D0_T_Matrix_
MueLu::DefaultLocalOrdinal LocalOrdinal
Teuchos::RCP< Teuchos::ParameterList > A22_RAP_reuse_data_
int BCedges_
Vectors for BCs.
void solveH(const MultiVector &RHS, MultiVector &X) const
apply solve to 1-1 block only
Teuchos::RCP< const Map > D0origDomainMap_
Teuchos::RCP< MultiVector > D0xSubComm_
Teuchos::RCP< Matrix > Ms_Matrix_
Teuchos::RCP< MultiVector > NullspaceH_
Nullspace for (1,1) problem.
RefMaxwell(const Teuchos::RCP< Matrix > &SM_Matrix, const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &M0inv_Matrix, const Teuchos::RCP< Matrix > &M1_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List, bool ComputePrec=true)
Teuchos::RCP< Matrix > D0_Matrix_
Teuchos::RCP< SmootherPrototype > PostSmootherData_
Teuchos::RCP< Matrix > SM_Matrix_
Various matrices.
Teuchos::RCP< MultiVector > D0res_
RefMaxwell(const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &M1_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List)
Teuchos::ParameterList smootherList_
void setFineLevelSmoother()
Set the fine level smoother.
Teuchos::RCP< Matrix > A22_
Teuchos::RCP< MultiVector > D0resTmp_
Teuchos::RCP< MultiVector > D0TR11Tmp_
Namespace for MueLu classes and methods.
Teuchos::RCP< Matrix > Addon_Matrix_
void setParameters(Teuchos::ParameterList &list)
Set parameters.
void applyInverseAdditive(const MultiVector &RHS, MultiVector &X) const
apply additive algorithm for 2x2 solve
void buildProlongator()
Setup the prolongator for the (1,1)-block.
MueLu::DefaultNode Node
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
Teuchos::RCP< MultiVector > P11res_
Temporary memory.
Static class that holds the complete list of valid MueLu parameters.
Teuchos::RCP< const Import > D0origImporter_
void allocateMemory(int numVectors) const
allocate multivectors for solve
Teuchos::RCP< MultiVector > P11resSubComm_
Teuchos::RCP< Matrix > M0inv_Matrix_
Teuchos::RCP< MultiVector > residual_
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
Teuchos::RCP< MultiVector > D0x_
Teuchos::RCP< MultiVector > P11resTmp_
void dumpCoords(const RealValuedMultiVector &X, std::string name) const
dump out real-valued multivector
MueLu::DefaultScalar Scalar
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Teuchos::RCP< Teuchos::ParameterList > AH_AP_reuse_data_
Teuchos::ParameterList precList11_
RefMaxwell(const Teuchos::RCP< Matrix > &SM_Matrix, const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &M1_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List, bool ComputePrec)
Teuchos::RCP< MultiVector > D0resSubComm_
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
Teuchos::RCP< Matrix > P11_
Teuchos::RCP< Matrix > A_nodal_Matrix_
Teuchos::RCP< Hierarchy > Hierarchy22_
Teuchos::RCP< Teuchos::ParameterList > AH_RAP_reuse_data_
Teuchos::RCP< Matrix > M1_Matrix_
Teuchos::RCP< RealValuedMultiVector > CoordsH_
Teuchos::ArrayRCP< bool > BCrows_
Teuchos::ArrayRCP< bool > BCcols_
Xpetra::MultiVector< coordinateType, LO, GO, NO > RealValuedMultiVector
Teuchos::RCP< SmootherBase > PreSmoother_
void compute(bool reuse=false)
Setup the preconditioner.
Teuchos::RCP< MultiVector > P11xSubComm_
Teuchos::ArrayRCP< bool > BCdomain_
bool disable_addon_
Some options.
void initialize(const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &Ms_Matrix, const Teuchos::RCP< Matrix > &M0inv_Matrix, const Teuchos::RCP< Matrix > &M1_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List)
virtual ~RefMaxwell()
Destructor.
Teuchos::ParameterList parameterList_
Parameter lists.
Teuchos::RCP< RealValuedMultiVector > Coords_
Coordinates.
RefMaxwell(const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &M0inv_Matrix, const Teuchos::RCP< Matrix > &M1_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List)
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
Teuchos::RCP< SmootherPrototype > PreSmootherData_
Teuchos::RCP< MultiVector > P11x_
Teuchos::RCP< Matrix > R11_
void residual(const MultiVector &X, const MultiVector &B, MultiVector &R) const
Compute a residual R = B - (*this) * X.
const Teuchos::RCP< Matrix > & getJacobian() const
Returns Jacobian matrix SM.
void dump(const Matrix &A, std::string name) const
dump out matrix
Teuchos::RCP< Matrix > AH_
void formCoarseMatrix()
Compute P11^{T}*A*P11 efficiently.
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
Teuchos::ScalarTraits< Scalar >::coordinateType coordinateType
Teuchos::RCP< Hierarchy > HierarchyH_
Two hierarchies: one for the coarse (1,1)-block, another for the (2,2)-block.
RefMaxwell(Teuchos::RCP< Hierarchy > HH, Teuchos::RCP< Hierarchy > H22)
Constructor with Hierarchies.
void solve22(const MultiVector &RHS, MultiVector &X) const
apply solve to 2-2 block only
Teuchos::RCP< Teuchos::ParameterList > A22_AP_reuse_data_
Teuchos::ParameterList precList22_
Teuchos::RCP< MultiVector > Nullspace_
Nullspace.
Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
RefMaxwell(const Teuchos::RCP< Matrix > &SM_Matrix, Teuchos::ParameterList &List, bool ComputePrec=true)
Teuchos::RCP< const Import > ImporterH_
Importer to coarse (1,1) hierarchy.
RefMaxwell(const Teuchos::RCP< Matrix > &SM_Matrix, const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &Ms_Matrix, const Teuchos::RCP< Matrix > &M0inv_Matrix, const Teuchos::RCP< Matrix > &M1_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List, bool ComputePrec=true)