MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_ShiftedLaplacian_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 // Jeremie Gaidamour (jngaida@sandia.gov)
40 // Jonathan Hu (jhu@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_SHIFTEDLAPLACIAN_DECL_HPP
47 #define MUELU_SHIFTEDLAPLACIAN_DECL_HPP
48 
49 // Xpetra
50 #include <Xpetra_Matrix_fwd.hpp>
53 #include <Xpetra_TpetraMultiVector.hpp>
54 
55 // MueLu
56 #include "MueLu.hpp"
57 #include "MueLu_ConfigDefs.hpp"
58 
59 #if defined(HAVE_MUELU_IFPACK2) and defined(HAVE_MUELU_TPETRA)
60 
61 #include <MueLu_BaseClass.hpp>
68 #include <MueLu_Hierarchy_fwd.hpp>
70 #include <MueLu_PFactory_fwd.hpp>
71 #include <MueLu_PgPFactory_fwd.hpp>
72 #include <MueLu_RAPFactory_fwd.hpp>
74 #include <MueLu_SaPFactory_fwd.hpp>
76 #include <MueLu_ShiftedLaplacianOperator.hpp>
82 #include <MueLu_Utilities_fwd.hpp>
83 
84 // Belos
85 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
86 #include <BelosConfigDefs.hpp>
87 #include <BelosLinearProblem.hpp>
88 #include <BelosSolverFactory.hpp>
89 #include <BelosTpetraAdapter.hpp>
90 #endif
91 
92 namespace MueLu {
93 
103  template <class Scalar = Xpetra::Matrix<>::scalar_type,
104  class LocalOrdinal = typename Xpetra::Matrix<Scalar>::local_ordinal_type,
105  class GlobalOrdinal = typename Xpetra::Matrix<Scalar, LocalOrdinal>::global_ordinal_type,
107  class ShiftedLaplacian : public BaseClass {
108 #undef MUELU_SHIFTEDLAPLACIAN_SHORT
109 #include "MueLu_UseShortNames.hpp"
110 
111  typedef Tpetra::Vector<SC,LO,GO,NO> TVEC;
112  typedef Tpetra::MultiVector<SC,LO,GO,NO> TMV;
113  typedef Tpetra::Operator<SC,LO,GO,NO> OP;
114 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
115  typedef Belos::LinearProblem<SC,TMV,OP> LinearProblem;
116  typedef Belos::SolverManager<SC,TMV,OP> SolverManager;
117  typedef Belos::SolverFactory<SC,TMV,OP> SolverFactory;
118 #endif
119 
120  public:
121 
122  /*
123  FIXME 26-June-2015 JJH: This contructor is setting numerous defaults. However, they don't match the defaults
124  FIXME int the method setParameters(). There also isn't any parameter validation that I can see.
125  */
126 
129  numPDEs_(1),
130  Smoother_("schwarz"),
131  Aggregation_("uncoupled"),
132  Nullspace_("constant"),
133  numLevels_(5),
134  coarseGridSize_(100),
135  omega_(2.0*M_PI),
136  iters_(500),
137  blksize_(1),
138  tol_(1.0e-4),
139  nsweeps_(5),
140  ncycles_(1),
141  cycles_(8),
142  subiters_(10),
143  option_(1),
144  nproblems_(0),
145  solverType_(1),
146  restart_size_(100),
147  recycle_size_(25),
148  smoother_sweeps_(4),
149  smoother_damping_((SC)1.0),
150  krylov_type_(1),
153  ilu_leveloffill_(5.0),
154  ilu_abs_thresh_(0.0),
155  ilu_rel_thresh_(1.0),
157  ilu_drop_tol_(0.01),
158  ilu_fill_tol_(0.01),
159  ilu_relax_val_(1.0),
160  ilu_rowperm_("LargeDiag"),
161  ilu_colperm_("COLAMD"),
162  ilu_drop_rule_("DROP_BASIC"),
163  ilu_normtype_("INF_NORM"),
164  ilu_milutype_("SILU"),
165  schwarz_overlap_(0),
166  schwarz_usereorder_(true),
167  schwarz_combinemode_(Tpetra::ADD),
168  schwarz_ordermethod_("rcm"),
169  GridTransfersExist_(false),
170  isSymmetric_(true)
171  { }
172 
173  // Destructor
174  virtual ~ShiftedLaplacian();
175 
176  // Parameters
178 
179  // Set matrices
180  void setProblemMatrix(RCP<Matrix>& A);
181  void setProblemMatrix(RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> >& TpetraA);
183  void setPreconditioningMatrix(RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> >& TpetraP);
184  void setstiff(RCP<Matrix>& K);
185  void setstiff(RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> >& TpetraK);
186  void setmass(RCP<Matrix>& M);
187  void setmass(RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> >& TpetraM);
188  void setcoords(RCP<MultiVector>& Coords);
189  void setNullSpace(RCP<MultiVector> NullSpace);
190  void setLevelShifts(std::vector<Scalar> levelshifts);
191 
192  // initialize: set parameters and factories, construct
193  // prolongation and restriction matrices
194  void initialize();
195  // setupFastRAP: setup hierarchy with
196  // prolongators of the stiffness matrix
197  // constant complex shifts
198  void setupFastRAP();
199  // setupSlowRAP: setup hierarchy with
200  // prolongators of the stiffness matrix
201  // variable complex shifts
202  void setupSlowRAP();
203  // setupNormalRAP: setup hierarchy with
204  // prolongators of the preconditioning matrix
205  void setupNormalRAP();
206  // setupSolver: initialize Belos solver
207  void setupSolver();
208  // resetLinearProblem: for multiple frequencies;
209  // reset the Belos operator if the frequency changes
210  void resetLinearProblem();
211 
212 
213  // Solve phase
214  int solve(const RCP<TMV> B, RCP<TMV>& X);
215  void multigrid_apply(const RCP<MultiVector> B,
216  RCP<MultiVector>& X);
217  void multigrid_apply(const RCP<Tpetra::MultiVector<SC,LO,GO,NO> > B,
218  RCP<Tpetra::MultiVector<SC,LO,GO,NO> >& X);
219  int GetIterations();
221 
223 
224  private:
225 
226  // Problem options
227  // numPDEs_ -> number of DOFs at each node
228  int numPDEs_;
229 
230  // Multigrid options
231  // numLevels_ -> number of Multigrid levels
232  // coarseGridSize_ -> size of coarsest grid (if current level has less DOFs, stop coarsening)
235 
236  // Shifted Laplacian/Helmholtz parameters
237  double omega_;
238  std::vector<SC> levelshifts_;
239 
240  // Krylov solver inputs
241  // iters -> max number of iterations
242  // tol -> residual tolerance
244  double tol_;
248 
249  // Smoother parameters
260  Tpetra::CombineMode schwarz_combinemode_;
261  std::string schwarz_ordermethod_;
262 
263  // flags for setup
266 
267  // Xpetra matrices
268  // K_ -> stiffness matrix
269  // M_ -> mass matrix
270  // A_ -> Problem matrix
271  // P_ -> Preconditioning matrix
274 
275  // Multigrid Hierarchy
277 
278  // Factories and prototypes
293  std::string precType_;
295 
296  // Operator and Preconditioner
299 
300 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
301  // Belos Linear Problem and Solver
302  RCP<LinearProblem> LinearProblem_;
303  RCP<SolverManager> SolverManager_;
304  RCP<SolverFactory> SolverFactory_;
305  RCP<Teuchos::ParameterList> BelosList_;
306 #endif
307 
308  };
309 
310 }
311 
312 #define MUELU_SHIFTEDLAPLACIAN_SHORT
313 
314 #endif //if defined(HAVE_MUELU_IFPACK2) and defined(HAVE_MUELU_TPETRA)
315 
316 #endif // MUELU_SHIFTEDLAPLACIAN_DECL_HPP
RCP< SmootherFactory > coarsestSmooFact_
void multigrid_apply(const RCP< MultiVector > B, RCP< MultiVector > &X)
RCP< CoarseMapFactory > CoarseMapfact_
RCP< CoalesceDropFactory > Dropfact_
void setParameters(Teuchos::RCP< Teuchos::ParameterList > paramList)
Teuchos::ScalarTraits< Scalar >::magnitudeType GetResidual()
Tpetra::MultiVector< SC, LO, GO, NO > TMV
RCP< SmootherPrototype > smooProto_
void setLevelShifts(std::vector< Scalar > levelshifts)
RCP< MueLu::ShiftedLaplacianOperator< SC, LO, GO, NO > > MueLuOp_
Teuchos::ParameterList coarsestSmooList_
void setcoords(RCP< MultiVector > &Coords)
RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > TpetraA_
RCP< UncoupledAggregationFactory > UCaggfact_
Tpetra::Operator< SC, LO, GO, NO > OP
int solve(const RCP< TMV > B, RCP< TMV > &X)
Base class for MueLu classes.
RCP< SmootherPrototype > coarsestSmooProto_
RCP< CoupledAggregationFactory > Aggfact_
Tpetra::Vector< SC, LO, GO, NO > TVEC
void setProblemMatrix(RCP< Matrix > &A)
Scalar SC
RCP< TentativePFactory > TentPfact_
Shifted Laplacian Helmholtz solver.
void setNullSpace(RCP< MultiVector > NullSpace)
void setPreconditioningMatrix(RCP< Matrix > &P)