All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
parallelamgbackend.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef EWOMS_PARALLEL_AMG_BACKEND_HH
28 #define EWOMS_PARALLEL_AMG_BACKEND_HH
29 
30 #include "parallelbasebackend.hh"
31 #include "bicgstabsolver.hh"
32 #include "combinedcriterion.hh"
33 
34 #include <dune/istl/paamg/amg.hh>
35 #include <dune/istl/paamg/pinfo.hh>
36 #include <dune/istl/owneroverlapcopy.hh>
37 
38 #include <iostream>
39 
40 namespace Ewoms {
41 namespace Linear {
42 template <class TypeTag>
44 }
45 
46 namespace Properties {
47 NEW_TYPE_TAG(ParallelAmgLinearSolver, INHERITS_FROM(ParallelBaseLinearSolver));
48 
49 NEW_PROP_TAG(AmgCoarsenTarget);
50 NEW_PROP_TAG(LinearSolverMaxError);
51 
54 SET_INT_PROP(ParallelAmgLinearSolver, AmgCoarsenTarget, 5000);
55 
56 SET_SCALAR_PROP(ParallelAmgLinearSolver, LinearSolverMaxError, 1e7);
57 
58 SET_TYPE_PROP(ParallelAmgLinearSolver, LinearSolverBackend,
60 } // namespace Properties
61 
62 namespace Linear {
69 template <class TypeTag>
70 class ParallelAmgBackend : public ParallelBaseBackend<TypeTag>
71 {
72  typedef ParallelBaseBackend<TypeTag> ParentType;
73 
74  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
75  typedef typename GET_PROP_TYPE(TypeTag, LinearSolverScalar) LinearSolverScalar;
76  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
77  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
78  typedef typename GET_PROP_TYPE(TypeTag, Overlap) Overlap;
79 
80  typedef typename ParentType::ParallelOperator ParallelOperator;
81  typedef typename ParentType::OverlappingVector OverlappingVector;
82  typedef typename ParentType::ParallelPreconditioner ParallelPreconditioner;
83  typedef typename ParentType::ParallelScalarProduct ParallelScalarProduct;
84 
85  static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq);
86  typedef Dune::FieldVector<LinearSolverScalar, numEq> VectorBlock;
87  typedef Dune::FieldMatrix<LinearSolverScalar, numEq, numEq> MatrixBlock;
88 
89  typedef Dune::BCRSMatrix<MatrixBlock> Matrix;
90  typedef Dune::BlockVector<VectorBlock> Vector;
91 
92  // define the smoother used for the AMG and specify its
93  // arguments
94  typedef Dune::SeqSOR<Matrix, Vector, Vector> SequentialSmoother;
95 // typedef Dune::SeqSSOR<Matrix,Vector,Vector> SequentialSmoother;
96 // typedef Dune::SeqJac<Matrix,Vector,Vector> SequentialSmoother;
97 // typedef Dune::SeqILU0<Matrix,Vector,Vector> SequentialSmoother;
98 // typedef Dune::SeqILUn<Matrix,Vector,Vector> SequentialSmoother;
99 
100 #if HAVE_MPI
101  typedef Dune::OwnerOverlapCopyCommunication<Ewoms::Linear::Index>
102  OwnerOverlapCopyCommunication;
103  typedef Dune::OverlappingSchwarzOperator<Matrix,
104  Vector,
105  Vector,
106  OwnerOverlapCopyCommunication> FineOperator;
107  typedef Dune::OverlappingSchwarzScalarProduct<Vector,
108  OwnerOverlapCopyCommunication> FineScalarProduct;
109  typedef Dune::BlockPreconditioner<Vector,
110  Vector,
111  OwnerOverlapCopyCommunication,
112  SequentialSmoother> ParallelSmoother;
113  typedef Dune::Amg::AMG<FineOperator,
114  Vector,
115  ParallelSmoother,
116  OwnerOverlapCopyCommunication> AMG;
117 #else
118  typedef Dune::MatrixAdapter<Matrix, Vector, Vector> FineOperator;
119  typedef Dune::SeqScalarProduct<Vector> FineScalarProduct;
120  typedef SequentialSmoother ParallelSmoother;
121  typedef Dune::Amg::AMG<FineOperator, Vector, ParallelSmoother> AMG;
122 #endif
123 
124  typedef BiCGStabSolver<ParallelOperator,
125  OverlappingVector,
126  AMG> RawLinearSolver;
127 
128 public:
129  ParallelAmgBackend(const Simulator& simulator)
130  : ParentType(simulator)
131  { }
132 
133  static void registerParameters()
134  {
136 
137  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverMaxError,
138  "The maximum residual error which the linear solver tolerates"
139  " without giving up");
140  EWOMS_REGISTER_PARAM(TypeTag, int, AmgCoarsenTarget,
141  "The coarsening target for the agglomerations of "
142  "the AMG preconditioner");
143  }
144 
145 protected:
146  friend ParentType;
147 
148  std::shared_ptr<AMG> preparePreconditioner_()
149  {
150 #if HAVE_MPI
151  // create and initialize DUNE's OwnerOverlapCopyCommunication
152  // using the domestic overlap
153  istlComm_ = std::make_shared<OwnerOverlapCopyCommunication>(MPI_COMM_WORLD);
154  setupAmgIndexSet_(this->overlappingMatrix_->overlap(), istlComm_->indexSet());
155  istlComm_->remoteIndices().template rebuild<false>();
156 #endif
157 
158  // create the parallel scalar product and the parallel operator
159 #if HAVE_MPI
160  fineOperator_ = std::make_shared<FineOperator>(*this->overlappingMatrix_, *istlComm_);
161 #else
162  fineOperator_ = std::make_shared<FineOperator>(*this->overlappingMatrix_);
163 #endif
164 
165  setupAmg_();
166 
167  return amg_;
168  }
169 
170  void cleanupPreconditioner_()
171  { /* nothing to do */ }
172 
173  std::shared_ptr<RawLinearSolver> prepareSolver_(ParallelOperator& parOperator,
174  ParallelScalarProduct& parScalarProduct,
175  AMG& parPreCond)
176  {
177  const auto& gridView = this->simulator_.gridView();
178  typedef CombinedCriterion<OverlappingVector, decltype(gridView.comm())> CCC;
179 
180  Scalar linearSolverTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance);
181  Scalar linearSolverAbsTolerance = this->simulator_.model().newtonMethod().tolerance() / 10.0;
182 
183  convCrit_.reset(new CCC(gridView.comm(),
184  /*residualReductionTolerance=*/linearSolverTolerance,
185  /*absoluteResidualTolerance=*/linearSolverAbsTolerance,
186  EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverMaxError)));
187 
188  auto bicgstabSolver =
189  std::make_shared<RawLinearSolver>(parPreCond, *convCrit_, parScalarProduct);
190 
191  int verbosity = 0;
192  if (parOperator.overlap().myRank() == 0)
193  verbosity = EWOMS_GET_PARAM(TypeTag, int, LinearSolverVerbosity);
194  bicgstabSolver->setVerbosity(verbosity);
195  bicgstabSolver->setMaxIterations(EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIterations));
196  bicgstabSolver->setLinearOperator(&parOperator);
197  bicgstabSolver->setRhs(this->overlappingb_);
198 
199  return bicgstabSolver;
200  }
201 
202  bool runSolver_(std::shared_ptr<RawLinearSolver> solver)
203  { return solver->apply(*this->overlappingx_); }
204 
205  void cleanupSolver_()
206  { /* nothing to do */ }
207 
208 #if HAVE_MPI
209  template <class ParallelIndexSet>
210  void setupAmgIndexSet_(const Overlap& overlap, ParallelIndexSet& istlIndices)
211  {
212  typedef Dune::OwnerOverlapCopyAttributeSet GridAttributes;
213  typedef Dune::OwnerOverlapCopyAttributeSet::AttributeSet GridAttributeSet;
214 
215  // create DUNE's ParallelIndexSet from a domestic overlap
216  istlIndices.beginResize();
217  for (Index curIdx = 0; static_cast<size_t>(curIdx) < overlap.numDomestic(); ++curIdx) {
218  GridAttributeSet gridFlag =
219  overlap.iAmMasterOf(curIdx)
220  ? GridAttributes::owner
221  : GridAttributes::copy;
222 
223  // an index is used by other processes if it is in the
224  // domestic or in the foreign overlap.
225  bool isShared = overlap.isInOverlap(curIdx);
226 
227  assert(curIdx == overlap.globalToDomestic(overlap.domesticToGlobal(curIdx)));
228  istlIndices.add(/*globalIdx=*/overlap.domesticToGlobal(curIdx),
229  Dune::ParallelLocalIndex<GridAttributeSet>(static_cast<size_t>(curIdx),
230  gridFlag,
231  isShared));
232  }
233  istlIndices.endResize();
234  }
235 #endif
236 
237  void setupAmg_()
238  {
239  if (amg_)
240  amg_.reset();
241 
242  int verbosity = 0;
243  if (this->simulator_.gridManager().gridView().comm().rank() == 0)
244  verbosity = EWOMS_GET_PARAM(TypeTag, int, LinearSolverVerbosity);
245 
246  typedef typename Dune::Amg::SmootherTraits<ParallelSmoother>::Arguments SmootherArgs;
247 
248  SmootherArgs smootherArgs;
249  smootherArgs.iterations = 1;
250  smootherArgs.relaxationFactor = 1.0;
251 
252  // specify the coarsen criterion:
253  //
254  // typedef
255  // Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Matrix,
256  // Dune::Amg::FirstDiagonal>>
257  typedef Dune::Amg::
258  CoarsenCriterion<Dune::Amg::SymmetricCriterion<Matrix, Dune::Amg::FrobeniusNorm> >
259  CoarsenCriterion;
260  int coarsenTarget = EWOMS_GET_PARAM(TypeTag, int, AmgCoarsenTarget);
261  CoarsenCriterion coarsenCriterion(/*maxLevel=*/15, coarsenTarget);
262  coarsenCriterion.setDefaultValuesAnisotropic(GridView::dimension,
263  /*aggregateSizePerDim=*/3);
264  if (verbosity > 0)
265  coarsenCriterion.setDebugLevel(1);
266  else
267  coarsenCriterion.setDebugLevel(0); // make the AMG shut up
268 
269  // reduce the minium coarsen rate (default is 1.2)
270  coarsenCriterion.setMinCoarsenRate(1.05);
271  // coarsenCriterion.setAccumulate(Dune::Amg::noAccu);
272  coarsenCriterion.setAccumulate(Dune::Amg::atOnceAccu);
273  coarsenCriterion.setSkipIsolated(false);
274 
275 // instantiate the AMG preconditioner
276 #if HAVE_MPI
277  amg_ = std::make_shared<AMG>(*fineOperator_, coarsenCriterion, smootherArgs, *istlComm_);
278 #else
279  amg_ = std::make_shared<AMG>(*fineOperator_, coarsenCriterion, smootherArgs);
280 #endif
281  }
282 
283  std::unique_ptr<ConvergenceCriterion<OverlappingVector> > convCrit_;
284 
285  std::shared_ptr<FineOperator> fineOperator_;
286  std::shared_ptr<AMG> amg_;
287 
288 #if HAVE_MPI
289  std::shared_ptr<OwnerOverlapCopyCommunication> istlComm_;
290 #endif
291 };
292 
293 } // namespace Linear
294 } // namespace Ewoms
295 
296 #endif
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
Provides a linear solver backend using the parallel algebraic multi-grid (AMG) linear solver from DUN...
Definition: parallelamgbackend.hh:43
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
GridManager & gridManager()
Return a reference to the grid manager of simulation.
Definition: simulator.hh:171
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
static void registerParameters()
Register all run-time parameters for the linear solver.
Definition: parallelbasebackend.hh:187
Implements a preconditioned stabilized BiCG linear solver.
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition: simulator.hh:183
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394
Provides the common code which is required by most linear solvers.