MueLu  Version of the Day
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>
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>
67 #include <MueLu_Hierarchy_fwd.hpp>
69 #include <MueLu_PFactory_fwd.hpp>
70 #include <MueLu_PgPFactory_fwd.hpp>
71 #include <MueLu_RAPFactory_fwd.hpp>
73 #include <MueLu_SaPFactory_fwd.hpp>
75 #include <MueLu_ShiftedLaplacianOperator.hpp>
81 #include <MueLu_Utilities_fwd.hpp>
82 
83 // Belos
84 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
85 #include <BelosConfigDefs.hpp>
86 #include <BelosLinearProblem.hpp>
87 #include <BelosSolverFactory.hpp>
88 #include <BelosTpetraAdapter.hpp>
89 #endif
90 
91 namespace MueLu {
92 
102  template <class Scalar = Xpetra::Matrix<>::scalar_type,
103  class LocalOrdinal = typename Xpetra::Matrix<Scalar>::local_ordinal_type,
104  class GlobalOrdinal = typename Xpetra::Matrix<Scalar, LocalOrdinal>::global_ordinal_type,
106  class ShiftedLaplacian : public BaseClass {
107 #undef MUELU_SHIFTEDLAPLACIAN_SHORT
108 #include "MueLu_UseShortNames.hpp"
109 
113 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
114  typedef Belos::LinearProblem<SC,TMV,OP> LinearProblem;
115  typedef Belos::SolverManager<SC,TMV,OP> SolverManager;
116  typedef Belos::SolverFactory<SC,TMV,OP> SolverFactory;
117 #endif
118 
119  public:
120 
121  /*
122  FIXME 26-June-2015 JJH: This contructor is setting numerous defaults. However, they don't match the defaults
123  FIXME int the method setParameters(). There also isn't any parameter validation that I can see.
124  */
125 
128  numPDEs_(1),
129  Smoother_("schwarz"),
130  Aggregation_("uncoupled"),
131  Nullspace_("constant"),
132  numLevels_(5),
133  coarseGridSize_(100),
134  omega_(2.0*M_PI),
135  iters_(500),
136  blksize_(1),
137  tol_(1.0e-4),
138  nsweeps_(5),
139  ncycles_(1),
140  cycles_(8),
141  subiters_(10),
142  option_(1),
143  nproblems_(0),
144  solverType_(1),
145  restart_size_(100),
146  recycle_size_(25),
147  smoother_sweeps_(4),
148  smoother_damping_((SC)1.0),
149  krylov_type_(1),
152  ilu_leveloffill_(5.0),
153  ilu_abs_thresh_(0.0),
154  ilu_rel_thresh_(1.0),
156  ilu_drop_tol_(0.01),
157  ilu_fill_tol_(0.01),
158  ilu_relax_val_(1.0),
159  ilu_rowperm_("LargeDiag"),
160  ilu_colperm_("COLAMD"),
161  ilu_drop_rule_("DROP_BASIC"),
162  ilu_normtype_("INF_NORM"),
163  ilu_milutype_("SILU"),
164  schwarz_overlap_(0),
165  schwarz_usereorder_(true),
167  schwarz_ordermethod_("rcm"),
168  GridTransfersExist_(false),
169  isSymmetric_(true)
170  { }
171 
172  // Destructor
173  virtual ~ShiftedLaplacian();
174 
175  // Parameters
176  void setParameters(Teuchos::RCP< Teuchos::ParameterList > paramList);
177 
178  // Set matrices
179  void setProblemMatrix(RCP<Matrix>& A);
180  void setProblemMatrix(RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> >& TpetraA);
181  void setPreconditioningMatrix(RCP<Matrix>& P);
183  void setstiff(RCP<Matrix>& K);
184  void setstiff(RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> >& TpetraK);
185  void setmass(RCP<Matrix>& M);
186  void setmass(RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> >& TpetraM);
187  void setcoords(RCP<MultiVector>& Coords);
188  void setNullSpace(RCP<MultiVector> NullSpace);
189  void setLevelShifts(std::vector<Scalar> levelshifts);
190 
191  // initialize: set parameters and factories, construct
192  // prolongation and restriction matrices
193  void initialize();
194  // setupFastRAP: setup hierarchy with
195  // prolongators of the stiffness matrix
196  // constant complex shifts
197  void setupFastRAP();
198  // setupSlowRAP: setup hierarchy with
199  // prolongators of the stiffness matrix
200  // variable complex shifts
201  void setupSlowRAP();
202  // setupNormalRAP: setup hierarchy with
203  // prolongators of the preconditioning matrix
204  void setupNormalRAP();
205  // setupSolver: initialize Belos solver
206  void setupSolver();
207  // resetLinearProblem: for multiple frequencies;
208  // reset the Belos operator if the frequency changes
209  void resetLinearProblem();
210 
211 
212  // Solve phase
213  int solve(const RCP<TMV> B, RCP<TMV>& X);
214  void multigrid_apply(const RCP<MultiVector> B,
215  RCP<MultiVector>& X);
218  int GetIterations();
219  double GetResidual();
220 
221  RCP<FactoryManager> Manager_;
222 
223  private:
224 
225  // Problem options
226  // numPDEs_ -> number of DOFs at each node
227  int numPDEs_;
228 
229  // Multigrid options
230  // numLevels_ -> number of Multigrid levels
231  // coarseGridSize_ -> size of coarsest grid (if current level has less DOFs, stop coarsening)
234 
235  // Shifted Laplacian/Helmholtz parameters
236  double omega_;
237  std::vector<SC> levelshifts_;
238 
239  // Krylov solver inputs
240  // iters -> max number of iterations
241  // tol -> residual tolerance
243  double tol_;
247 
248  // Smoother parameters
260  std::string schwarz_ordermethod_;
261 
262  // flags for setup
265 
266  // Xpetra matrices
267  // K_ -> stiffness matrix
268  // M_ -> mass matrix
269  // A_ -> Problem matrix
270  // P_ -> Preconditioning matrix
271  RCP<Matrix> K_, M_, A_, P_;
272  RCP<MultiVector> Coords_, NullSpace_;
273 
274  // Multigrid Hierarchy
275  RCP<Hierarchy> Hierarchy_;
276 
277  // Factories and prototypes
278  RCP<TentativePFactory> TentPfact_;
279  RCP<PFactory> Pfact_;
280  RCP<PgPFactory> PgPfact_;
281  RCP<TransPFactory> TransPfact_;
282  RCP<GenericRFactory> Rfact_;
283  RCP<RAPFactory> Acfact_;
284  RCP<RAPShiftFactory> Acshift_;
285  RCP<CoalesceDropFactory> Dropfact_;
286  RCP<CoupledAggregationFactory> Aggfact_;
287  RCP<UncoupledAggregationFactory> UCaggfact_;
288  RCP<SmootherPrototype> smooProto_, coarsestSmooProto_;
289  RCP<SmootherFactory> smooFact_, coarsestSmooFact_;
290  Teuchos::ParameterList coarsestSmooList_;
291  std::string precType_;
292  Teuchos::ParameterList precList_;
293 
294  // Operator and Preconditioner
295  RCP< MueLu::ShiftedLaplacianOperator<SC,LO,GO,NO> > MueLuOp_;
296  RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > TpetraA_;
297 
298 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
299  // Belos Linear Problem and Solver
300  RCP<LinearProblem> LinearProblem_;
301  RCP<SolverManager> SolverManager_;
302  RCP<SolverFactory> SolverFactory_;
303  RCP<Teuchos::ParameterList> BelosList_;
304 #endif
305 
306  };
307 
308 }
309 
310 #define MUELU_SHIFTEDLAPLACIAN_SHORT
311 
312 #endif //if defined(HAVE_MUELU_IFPACK2) and defined(HAVE_MUELU_TPETRA)
313 
314 #endif // MUELU_SHIFTEDLAPLACIAN_DECL_HPP
RCP< SmootherFactory > coarsestSmooFact_
void multigrid_apply(const RCP< MultiVector > B, RCP< MultiVector > &X)
LocalOrdinal local_ordinal_type
Namespace for MueLu classes and methods.
RCP< CoalesceDropFactory > Dropfact_
void setParameters(Teuchos::RCP< Teuchos::ParameterList > paramList)
Tpetra::MultiVector< SC, LO, GO, NO > TMV
RCP< SmootherPrototype > smooProto_
void setLevelShifts(std::vector< Scalar > levelshifts)
ADD
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)
GlobalOrdinal global_ordinal_type