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"
52 
58 #include "MueLu_Utilities_fwd.hpp"
59 
60 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
67 #endif
68 
70 #include "MueLu_TrilinosSmoother.hpp"
71 #include "MueLu_Hierarchy.hpp"
72 
73 #include "Xpetra_Map_fwd.hpp"
74 #include "Xpetra_Matrix_fwd.hpp"
75 #include "Xpetra_MatrixFactory_fwd.hpp"
76 #include "Xpetra_MultiVectorFactory_fwd.hpp"
77 #include "Xpetra_VectorFactory_fwd.hpp"
78 #include "Xpetra_CrsMatrixWrap_fwd.hpp"
79 #include "Xpetra_ExportFactory_fwd.hpp"
80 
81 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || defined(HAVE_MUELU_INST_DOUBLE_INT_INT))
82 #include "Ifpack2_Preconditioner.hpp"
83 #include "Ifpack2_Hiptmair.hpp"
84 #endif
85 
86 namespace MueLu {
87 
112  template <class Scalar,
113  class LocalOrdinal,
114  class GlobalOrdinal,
115  class Node>
116  class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
117 
118 #undef MUELU_REFMAXWELL_SHORT
119 #include "MueLu_UseShortNames.hpp"
120 
121  public:
122 
123  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
124  typedef typename Xpetra::MultiVector<magnitudeType,LO,GO,NO> RealValuedMultiVector;
125 
128  HierarchyH_(Teuchos::null),
129  Hierarchy22_(Teuchos::null),
130  disable_addon_(true),
131  mode_("additive")
132  {
133  }
134 
136  RefMaxwell(Teuchos::RCP<Hierarchy> HH, Teuchos::RCP<Hierarchy> H22) :
137  HierarchyH_(HH),
138  Hierarchy22_(H22),
139  disable_addon_(false),
140  mode_("additive")
141  {
142  }
143 
155  RefMaxwell(const Teuchos::RCP<Matrix> & SM_Matrix,
156  const Teuchos::RCP<Matrix> & D0_Matrix,
157  const Teuchos::RCP<Matrix> & M0inv_Matrix,
158  const Teuchos::RCP<Matrix> & M1_Matrix,
159  const Teuchos::RCP<MultiVector> & Nullspace,
160  const Teuchos::RCP<RealValuedMultiVector> & Coords,
161  Teuchos::ParameterList& List,
162  bool ComputePrec = true)
163  {
164  initialize(D0_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
165 
166  resetMatrix(SM_Matrix);
167 
168  // compute preconditioner (optionally)
169  if(ComputePrec)
170  compute();
171  }
172 
182  RefMaxwell(const Teuchos::RCP<Matrix> & D0_Matrix,
183  const Teuchos::RCP<Matrix> & M0inv_Matrix,
184  const Teuchos::RCP<Matrix> & M1_Matrix,
185  const Teuchos::RCP<MultiVector> & Nullspace,
186  const Teuchos::RCP<RealValuedMultiVector> & Coords,
187  Teuchos::ParameterList& List) : SM_Matrix_(Teuchos::null)
188  {
189  initialize(D0_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
190  }
191 
202  RefMaxwell(const Teuchos::RCP<Matrix> & SM_Matrix,
203  const Teuchos::RCP<Matrix> & D0_Matrix,
204  const Teuchos::RCP<Matrix> & M1_Matrix,
205  const Teuchos::RCP<MultiVector> & Nullspace,
206  const Teuchos::RCP<RealValuedMultiVector> & Coords,
207  Teuchos::ParameterList& List,
208  bool ComputePrec = true)
209  {
210  initialize(D0_Matrix,Teuchos::null,M1_Matrix,Nullspace,Coords,List);
211 
212  resetMatrix(SM_Matrix);
213 
214  // compute preconditioner (optionally)
215  if(ComputePrec)
216  compute();
217  }
218 
227  RefMaxwell(const Teuchos::RCP<Matrix> & D0_Matrix,
228  const Teuchos::RCP<Matrix> & M1_Matrix,
229  const Teuchos::RCP<MultiVector> & Nullspace,
230  const Teuchos::RCP<RealValuedMultiVector> & Coords,
231  Teuchos::ParameterList& List) : SM_Matrix_(Teuchos::null)
232  {
233  initialize(D0_Matrix,Teuchos::null,M1_Matrix,Nullspace,Coords,List);
234  }
235 
242  RefMaxwell(const Teuchos::RCP<Matrix> & SM_Matrix,
243  Teuchos::ParameterList& List,
244  bool ComputePrec = true)
245  {
246 
247  RCP<MultiVector> Nullspace = List.get<RCP<MultiVector> >("Nullspace", Teuchos::null);
248  RCP<RealValuedMultiVector> Coords = List.get<RCP<RealValuedMultiVector> >("Coordinates", Teuchos::null);
249  RCP<Matrix> D0_Matrix = List.get<RCP<Matrix> >("D0");
250  RCP<Matrix> M1_Matrix = List.get<RCP<Matrix> >("M1");
251  RCP<Matrix> M0inv_Matrix = List.get<RCP<Matrix> >("M0inv", Teuchos::null);
252 
253  initialize(D0_Matrix,M0inv_Matrix,M1_Matrix,Nullspace,Coords,List);
254 
255  if (SM_Matrix != Teuchos::null)
256  resetMatrix(SM_Matrix);
257 
258  // compute preconditioner (optionally)
259  if(ComputePrec)
260  compute();
261  }
262 
264  virtual ~RefMaxwell() {}
265 
267  Teuchos::RCP<const Map> getDomainMap() const;
268 
270  Teuchos::RCP<const Map> getRangeMap() const;
271 
273  const Teuchos::RCP<Matrix> & getJacobian() const {
274  return SM_Matrix_;
275  }
276 
278  void setParameters(Teuchos::ParameterList& list);
279 
281  void compute();
282 
284  void buildProlongator();
285 
287  void formCoarseMatrix();
288 
290  void resetMatrix(Teuchos::RCP<Matrix> SM_Matrix_new);
291 
293  void applyInverseAdditive(const MultiVector& RHS, MultiVector& X) const;
294 
296  void applyInverse121(const MultiVector& RHS, MultiVector& X) const;
297 
299  void applyInverse212(const MultiVector& RHS, MultiVector& X) const;
300 
302  void applyInverse11only(const MultiVector& RHS, MultiVector& X) const;
303 
307  void apply (const MultiVector& X, MultiVector& Y,
308  Teuchos::ETransp mode = Teuchos::NO_TRANS,
309  Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
310  Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const;
311 
313  bool hasTransposeApply() const;
314 
315  template <class NewNode>
316  Teuchos::RCP< RefMaxwell<Scalar, LocalOrdinal, GlobalOrdinal, NewNode> >
317  clone (const RCP<NewNode>& new_node) const {
319  (HierarchyH_->template clone<NewNode> (new_node),
320  Hierarchy22_->template clone<NewNode> (new_node)));
321  }
322 
323  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_HIGH) const;
324 
325  private:
326 
336  void initialize(const Teuchos::RCP<Matrix> & D0_Matrix,
337  const Teuchos::RCP<Matrix> & M0inv_Matrix,
338  const Teuchos::RCP<Matrix> & M1_Matrix,
339  const Teuchos::RCP<MultiVector> & Nullspace,
340  const Teuchos::RCP<RealValuedMultiVector> & Coords,
341  Teuchos::ParameterList& List);
342 
344  Teuchos::RCP<Hierarchy> HierarchyH_, Hierarchy22_;
345  Teuchos::RCP<SmootherBase> PreSmoother_, PostSmoother_;
346 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || defined(HAVE_MUELU_INST_DOUBLE_INT_INT))
347  Teuchos::RCP<Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> > hiptmairPreSmoother_, hiptmairPostSmoother_;
348 #endif
352  Teuchos::RCP<Matrix> A_nodal_Matrix_, P11_, R11_, AH_, A22_;
354 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
355  Kokkos::View<const bool*, typename Node::device_type> BCrowsKokkos_;
356  Kokkos::View<const bool*, typename Node::device_type> BCcolsKokkos_;
357 #endif
358  Teuchos::ArrayRCP<const bool> BCrows_;
359  Teuchos::ArrayRCP<const bool> BCcols_;
361  Teuchos::RCP<MultiVector> Nullspace_;
363  Teuchos::RCP<RealValuedMultiVector> Coords_, CoordsH_;
368  std::string mode_;
370  mutable Teuchos::RCP<MultiVector> P11res_, P11x_, D0res_, D0x_, residual_;
371  };
372 
373 } // namespace
374 
375 #define MUELU_REFMAXWELL_SHORT
376 #endif // MUELU_REFMAXWELL_DECL_HPP
Teuchos::RCP< SmootherBase > PostSmoother_
Teuchos::RCP< Matrix > D0_T_Matrix_
Teuchos::ArrayRCP< const bool > BCrows_
Vectors for BCs.
Teuchos::RCP< Matrix > Ms_Matrix_
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< 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_
Teuchos::ArrayRCP< const bool > BCcols_
Teuchos::RCP< Matrix > A22_
Namespace for MueLu classes and methods.
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.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
Teuchos::RCP< MultiVector > P11res_
Temporary memory.
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=true)
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell&#39;s equations in curl-curl form...
Xpetra::MultiVector< magnitudeType, LO, GO, NO > RealValuedMultiVector
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_
void applyInverse121(const MultiVector &RHS, MultiVector &X) const
apply 1-2-1 algorithm for 2x2 solve
Teuchos::ParameterList precList11_
void compute()
Setup the preconditioner.
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< Matrix > M1_Matrix_
Teuchos::RCP< RealValuedMultiVector > CoordsH_
void initialize(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)
Teuchos::RCP< SmootherBase > PreSmoother_
bool disable_addon_
Some options.
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.
void applyInverse212(const MultiVector &RHS, MultiVector &X) const
apply 2-1-2 algorithm for 2x2 solve
Teuchos::RCP< RefMaxwell< Scalar, LocalOrdinal, GlobalOrdinal, NewNode > > clone(const RCP< NewNode > &new_node) const
Teuchos::RCP< MultiVector > P11x_
Teuchos::RCP< Matrix > R11_
const Teuchos::RCP< Matrix > & getJacobian() const
Returns Jacobian matrix SM.
Teuchos::RCP< Matrix > AH_
void formCoarseMatrix()
Compute P11^{T}*A*P11 efficiently.
void applyInverse11only(const MultiVector &RHS, MultiVector &X) const
apply solve to 1-1 block only
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::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.
Teuchos::ParameterList precList22_
Teuchos::RCP< MultiVector > Nullspace_
Nullspace.
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new)
Reset system matrix.
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)