MueLu  Version of the Day
MueLu_Maxwell1_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_MAXWELL1_DECL_HPP
47 #define MUELU_MAXWELL1_DECL_HPP
48 
49 #include "MueLu_ConfigDefs.hpp"
50 #include "MueLu_BaseClass.hpp"
51 
53 #include "MueLu_SaPFactory_fwd.hpp"
54 
56 #include "MueLu_TrilinosSmoother.hpp"
57 #include "MueLu_Utilities_fwd.hpp"
58 #include "MueLu_Level_fwd.hpp"
59 #include "MueLu_Hierarchy_fwd.hpp"
60 #include "MueLu_RAPFactory_fwd.hpp"
61 #include "MueLu_PerfUtils_fwd.hpp"
62 #include "MueLu_SmootherBase.hpp"
63 
64 #if defined(HAVE_MUELU_KOKKOS_REFACTOR)
66 #endif
67 
68 #include "Xpetra_Map_fwd.hpp"
69 #include "Xpetra_Matrix_fwd.hpp"
70 #include "Xpetra_MatrixFactory_fwd.hpp"
71 #include "Xpetra_MultiVectorFactory_fwd.hpp"
72 #include "Xpetra_VectorFactory_fwd.hpp"
73 #include "Xpetra_CrsMatrixWrap_fwd.hpp"
74 
75 namespace MueLu {
76 
83  template <class Scalar,
84  class LocalOrdinal,
85  class GlobalOrdinal,
86  class Node>
87  class Maxwell1 : public VerboseObject, public Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
88 
89 #undef MUELU_MAXWELL1_SHORT
90 #include "MueLu_UseShortNames.hpp"
91 
92  public:
93 
94  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
95  typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType coordinateType;
96  typedef typename Xpetra::MultiVector<coordinateType,LO,GO,NO> RealValuedMultiVector;
97 
100  Hierarchy11_(Teuchos::null),
101  Hierarchy22_(Teuchos::null)
102  {
103  }
104 
106  Maxwell1(Teuchos::RCP<Hierarchy> H11, Teuchos::RCP<Hierarchy> H22) :
107  Hierarchy11_(H11),
108  Hierarchy22_(H22)
109  {
110  }
111 
121  Maxwell1(const Teuchos::RCP<Matrix> & SM_Matrix,
122  const Teuchos::RCP<Matrix> & D0_Matrix,
123  const Teuchos::RCP<MultiVector> & Nullspace,
124  const Teuchos::RCP<RealValuedMultiVector> & Coords,
125  Teuchos::ParameterList& List,
126  bool ComputePrec = true)
127  {
128  RCP<Matrix> Kn_Matrix;
129  initialize(D0_Matrix,Kn_Matrix,Nullspace,Coords,List);
130  resetMatrix(SM_Matrix,ComputePrec);
131  }
132 
142  Maxwell1(const Teuchos::RCP<Matrix> & SM_Matrix,
143  const Teuchos::RCP<Matrix> & D0_Matrix,
144  const Teuchos::RCP<Matrix> & Kn_Matrix,
145  const Teuchos::RCP<MultiVector> & Nullspace,
146  const Teuchos::RCP<RealValuedMultiVector> & Coords,
147  Teuchos::ParameterList& List,
148  bool ComputePrec = true)
149  {
150  initialize(D0_Matrix,Kn_Matrix,Nullspace,Coords,List);
151  resetMatrix(SM_Matrix,ComputePrec);
152  }
153 
154 
161  Maxwell1(const Teuchos::RCP<Matrix> & SM_Matrix,
162  Teuchos::ParameterList& List,
163  bool ComputePrec = true)
164  {
165 
166  RCP<MultiVector> Nullspace = List.get<RCP<MultiVector> >("Nullspace", Teuchos::null);
167  RCP<RealValuedMultiVector> Coords = List.get<RCP<RealValuedMultiVector> >("Coordinates", Teuchos::null);
168  RCP<Matrix> D0_Matrix = List.get<RCP<Matrix> >("D0");
169  RCP<Matrix> Kn_Matrix;
170  if (List.isType<RCP<Matrix> >("Kn"))
171  Kn_Matrix = List.get<RCP<Matrix> >("Kn");
172 
173  initialize(D0_Matrix,Kn_Matrix,Nullspace,Coords,List);
174 
175  if (SM_Matrix != Teuchos::null)
176  resetMatrix(SM_Matrix,ComputePrec);
177  }
178 
180  virtual ~Maxwell1() {}
181 
183  Teuchos::RCP<const Map> getDomainMap() const;
184 
186  Teuchos::RCP<const Map> getRangeMap() const;
187 
189  const Teuchos::RCP<Matrix> & getJacobian() const {
190  return SM_Matrix_;
191  }
192 
194  void setParameters(Teuchos::ParameterList& list);
195 
197  void compute(bool reuse=false);
198 
200  void resetMatrix(Teuchos::RCP<Matrix> SM_Matrix_new, bool ComputePrec=true);
201 
205  void apply (const MultiVector& X, MultiVector& Y,
206  Teuchos::ETransp mode = Teuchos::NO_TRANS,
207  Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
208  Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const;
209 
211  bool hasTransposeApply() const;
212 
213  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_HIGH) const;
214 
216  void residual(const MultiVector & X,
217  const MultiVector & B,
218  MultiVector & R) const {
219  using STS = Teuchos::ScalarTraits<Scalar>;
220  R.update(STS::one(),B,STS::zero());
221  this->apply (X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
222  }
223 
224 
225  private:
227  Teuchos::RCP<Matrix> generate_kn() const;
228 
229 
238  void initialize(const Teuchos::RCP<Matrix> & D0_Matrix,
239  const Teuchos::RCP<Matrix> & Kn_Matrix,
240  const Teuchos::RCP<MultiVector> & Nullspace,
241  const Teuchos::RCP<RealValuedMultiVector> & Coords,
242  Teuchos::ParameterList& List);
243 
244 
246  void applyInverseRefMaxwellAdditive(const MultiVector& RHS, MultiVector& X) const;
247 
249  void applyInverseStandard(const MultiVector& RHS, MultiVector& X) const;
250 
252  void allocateMemory(int numVectors) const;
253 
255  void dump(const Matrix& A, std::string name) const;
256 
258  void dump(const MultiVector& X, std::string name) const;
259 
261  void dumpCoords(const RealValuedMultiVector& X, std::string name) const;
262 
264  void dump(const Teuchos::ArrayRCP<bool>& v, std::string name) const;
265 
266 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
267  void dump(const Kokkos::View<bool*, typename Node::device_type>& v, std::string name) const;
269 #endif
270 
272  Teuchos::RCP<Teuchos::TimeMonitor> getTimer(std::string name, RCP<const Teuchos::Comm<int> > comm=Teuchos::null) const;
273 
275  mutable Teuchos::ParameterList parameterList_, precList11_, precList22_;
276 
278  Teuchos::RCP<Hierarchy> Hierarchy11_, Hierarchy22_;
279 
281  Teuchos::RCP<Matrix> SM_Matrix_, D0_Matrix_, Kn_Matrix_;
282 
284 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
285  Kokkos::View<bool*, typename Node::device_type> BCrowsKokkos_, BCcolsKokkos_, BCdomainKokkos_;
286 #endif
288  Teuchos::ArrayRCP<bool> BCrows_, BCcols_, BCdomain_;
290  Teuchos::RCP<MultiVector> Nullspace_;
292  Teuchos::RCP<RealValuedMultiVector> Coords_;
296 
300 
302  RCP<Matrix> P11_;
303  mutable Teuchos::RCP<MultiVector> residualFine_, residual11c_, residual22_, update11c_,update22_;
304 
305  };
306 
307 
308 
309 
310 } // namespace
311 
312 #define MUELU_MAXWELL1_SHORT
313 #endif // MUELU_MAXWELL1_DECL_HPP
Teuchos::RCP< MultiVector > Nullspace_
Nullspace.
MueLu::DefaultLocalOrdinal LocalOrdinal
Maxwell1(const Teuchos::RCP< Matrix > &SM_Matrix, const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List, bool ComputePrec=true)
Teuchos::RCP< Hierarchy > Hierarchy11_
Two hierarchies: one for the (1,1)-block, another for the (2,2)-block.
Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
void applyInverseStandard(const MultiVector &RHS, MultiVector &X) const
apply standard Maxwell1 cycle
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
void applyInverseRefMaxwellAdditive(const MultiVector &RHS, MultiVector &X) const
apply RefMaxwell additive 2x2 style cycle
bool useKokkos_
Some options.
Teuchos::ArrayRCP< bool > BCcols_
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
Teuchos::RCP< Matrix > D0_Matrix_
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)
Maxwell1(const Teuchos::RCP< Matrix > &SM_Matrix, 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, bool ComputePrec=true)
Teuchos::RCP< MultiVector > residual11c_
Teuchos::ParameterList precList22_
Namespace for MueLu classes and methods.
MueLu::DefaultNode Node
Teuchos::RCP< MultiVector > update22_
Verbose class for MueLu classes.
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell&#39;s equations in curl-curl form...
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
void setParameters(Teuchos::ParameterList &list)
Set parameters.
Teuchos::ArrayRCP< bool > BCdomain_
Teuchos::ParameterList precList11_
Teuchos::RCP< Hierarchy > Hierarchy22_
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Teuchos::RCP< Teuchos::TimeMonitor > getTimer(std::string name, RCP< const Teuchos::Comm< int > > comm=Teuchos::null) const
get a (synced) timer
Teuchos::ScalarTraits< Scalar >::coordinateType coordinateType
mode_type
Execution modes.
int BCedges_
Vectors for BCs.
Teuchos::RCP< Matrix > Kn_Matrix_
void dump(const Matrix &A, std::string name) const
dump out matrix
const Teuchos::RCP< Matrix > & getJacobian() const
Returns Jacobian matrix SM.
Maxwell1()
Constructor.
RCP< Matrix > P11_
Temporary memory (cached vectors for RefMaxwell-style)
Teuchos::ArrayRCP< bool > BCrows_
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
Maxwell1(Teuchos::RCP< Hierarchy > H11, Teuchos::RCP< Hierarchy > H22)
Constructor with Hierarchies.
void compute(bool reuse=false)
Setup the preconditioner.
Teuchos::RCP< Matrix > generate_kn() const
Generates the Kn matrix.
Xpetra::MultiVector< coordinateType, LO, GO, NO > RealValuedMultiVector
Teuchos::RCP< MultiVector > update11c_
Maxwell1(const Teuchos::RCP< Matrix > &SM_Matrix, Teuchos::ParameterList &List, bool ComputePrec=true)
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
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
virtual ~Maxwell1()
Destructor.
Teuchos::RCP< MultiVector > residual22_
void dumpCoords(const RealValuedMultiVector &X, std::string name) const
dump out real-valued multivector
Teuchos::RCP< RealValuedMultiVector > Coords_
Coordinates.
void residual(const MultiVector &X, const MultiVector &B, MultiVector &R) const
Compute a residual R = B - (*this) * X.
Teuchos::RCP< Matrix > SM_Matrix_
Various matrices.
Teuchos::RCP< MultiVector > residualFine_
Teuchos::ParameterList parameterList_
ParameterLists.
void allocateMemory(int numVectors) const
allocate multivectors for solve