parallelbasebackend.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_BASE_BACKEND_HH
28 #define EWOMS_PARALLEL_BASE_BACKEND_HH
29 
37 
41 
42 #include <dune/grid/io/file/vtk/vtkwriter.hh>
43 
44 #include <dune/common/fvector.hh>
45 
46 #include <sstream>
47 #include <memory>
48 #include <iostream>
49 
50 namespace Ewoms {
51 namespace Properties {
52 NEW_TYPE_TAG(ParallelBaseLinearSolver);
53 
54 // forward declaration of the required property tags
55 NEW_PROP_TAG(Simulator);
56 NEW_PROP_TAG(Scalar);
57 NEW_PROP_TAG(NumEq);
58 NEW_PROP_TAG(JacobianMatrix);
59 NEW_PROP_TAG(GlobalEqVector);
60 NEW_PROP_TAG(VertexMapper);
61 NEW_PROP_TAG(GridView);
62 
63 NEW_PROP_TAG(BorderListCreator);
64 NEW_PROP_TAG(Overlap);
65 NEW_PROP_TAG(OverlappingVector);
66 NEW_PROP_TAG(OverlappingMatrix);
67 NEW_PROP_TAG(OverlappingScalarProduct);
68 NEW_PROP_TAG(OverlappingLinearOperator);
69 
71 NEW_PROP_TAG(LinearSolverBackend);
72 
74 NEW_PROP_TAG(PreconditionerWrapper);
75 
76 
78 NEW_PROP_TAG(LinearSolverScalar);
79 
87 NEW_PROP_TAG(LinearSolverOverlapSize);
88 
92 NEW_PROP_TAG(LinearSolverTolerance);
93 
101 NEW_PROP_TAG(LinearSolverVerbosity);
102 
104 NEW_PROP_TAG(LinearSolverMaxIterations);
105 
107 NEW_PROP_TAG(PreconditionerOrder);
108 
110 NEW_PROP_TAG(PreconditionerRelaxation);
111 }} // namespace Properties, Ewoms
112 
113 namespace Ewoms {
114 namespace Linear {
141 template <class TypeTag>
143 {
144 protected:
145  typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) Implementation;
146 
147  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
148  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
149  typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) Matrix;
150  typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
151  typedef typename GET_PROP_TYPE(TypeTag, BorderListCreator) BorderListCreator;
152  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
153 
154  typedef typename GET_PROP_TYPE(TypeTag, Overlap) Overlap;
155  typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
156  typedef typename GET_PROP_TYPE(TypeTag, OverlappingMatrix) OverlappingMatrix;
157 
158  typedef typename GET_PROP_TYPE(TypeTag, PreconditionerWrapper) PreconditionerWrapper;
159  typedef typename PreconditionerWrapper::SequentialPreconditioner SequentialPreconditioner;
160 
161  typedef Ewoms::Linear::OverlappingPreconditioner<SequentialPreconditioner,
162  Overlap> ParallelPreconditioner;
163  typedef Ewoms::Linear::OverlappingScalarProduct<OverlappingVector,
164  Overlap> ParallelScalarProduct;
165  typedef Ewoms::Linear::OverlappingOperator<OverlappingMatrix,
166  OverlappingVector,
167  OverlappingVector> ParallelOperator;
168 
169  enum { dimWorld = GridView::dimensionworld };
170 
171 public:
172  ParallelBaseBackend(const Simulator& simulator)
173  : simulator_(simulator)
174  , gridSequenceNumber_( -1 )
175  {
176  overlappingMatrix_ = nullptr;
177  overlappingb_ = nullptr;
178  overlappingx_ = nullptr;
179  }
180 
182  { cleanup_(); }
183 
187  static void registerParameters()
188  {
189  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverTolerance,
190  "The maximum allowed error between of the linear solver");
191  EWOMS_REGISTER_PARAM(TypeTag, unsigned, LinearSolverOverlapSize,
192  "The size of the algebraic overlap for the linear solver");
193  EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIterations,
194  "The maximum number of iterations of the linear solver");
195  EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity,
196  "The verbosity level of the linear solver");
197 
198  PreconditionerWrapper::registerParameters();
199  }
200 
205  void eraseMatrix()
206  { cleanup_(); }
207 
208  void prepareMatrix(const Matrix& M)
209  {
210  // make sure that the overlapping matrix and block vectors
211  // have been created
212  prepare_(M);
213 
214  // copy the interior values of the non-overlapping linear system of
215  // equations to the overlapping one. On ther border, we add up
216  // the values of all processes (using the assignAdd() methods)
217  overlappingMatrix_->assignFromNative(M);
218 
219  asImp_().rescale_();
220 
221  // synchronize all entries from their master processes and add entries on the
222  // process border
223  overlappingMatrix_->syncAdd();
224  // the entries on the border have already been added in prepareRhs()
225  overlappingb_->sync();
226  }
227 
228  void prepareRhs(const Matrix& M, Vector& b)
229  {
230  // make sure that the overlapping matrix and block vectors
231  // have been created
232  prepare_(M);
233 
234  overlappingb_->assignAddBorder(b);
235 
236  // copy the result back to the non-overlapping vector. This is
237  // necessary here as assignAddBorder() might modify the
238  // residual vector for the border entities and we need the
239  // "globalized" residual in b...
240  overlappingb_->assignTo(b);
241  }
242 
248  bool solve(Vector& x)
249  {
250  (*overlappingx_) = 0.0;
251 
252  auto parPreCond = asImp_().preparePreconditioner_();
253 
254  auto cleanupPrecondFn =
255  [this]() -> void
256  { this->asImp_().cleanupPreconditioner_(); };
257 
258  GenericGuard<decltype(cleanupPrecondFn)> precondGuard(cleanupPrecondFn);
259 
260  // create the parallel scalar product and the parallel operator
261  ParallelScalarProduct parScalarProduct(overlappingMatrix_->overlap());
262  ParallelOperator parOperator(*overlappingMatrix_);
263 
264  // retrieve the linear solver
265  auto solver = asImp_().prepareSolver_(parOperator,
266  parScalarProduct,
267  *parPreCond);
268 
269  auto cleanupSolverFn =
270  [this]() -> void
271  { this->asImp_().cleanupSolver_(); };
272  GenericGuard<decltype(cleanupSolverFn)> solverGuard(cleanupSolverFn);
273 
274  // run the linear solver and have some fun
275  bool result = asImp_().runSolver_(solver);
276 
277  // copy the result back to the non-overlapping vector
278  overlappingx_->assignTo(x);
279 
280  // return the result of the solver
281  return result;
282  }
283 
284 protected:
285  Implementation& asImp_()
286  { return *static_cast<Implementation *>(this); }
287 
288  const Implementation& asImp_() const
289  { return *static_cast<const Implementation *>(this); }
290 
291  void prepare_(const Matrix& M)
292  {
293  // if grid has changed the sequence number has changed too
294  int curSeqNum = simulator_.gridManager().gridSequenceNumber();
295  if( gridSequenceNumber_ == curSeqNum && overlappingMatrix_)
296  // the grid has not changed since the overlappingMatrix_has been created, so
297  // there's noting to do
298  return;
299 
300  asImp_().cleanup_();
301  gridSequenceNumber_ = curSeqNum;
302 
303  BorderListCreator borderListCreator(simulator_.gridView(),
304  simulator_.model().dofMapper());
305 
306  // create the overlapping Jacobian matrix
307  unsigned overlapSize = EWOMS_GET_PARAM(TypeTag, unsigned, LinearSolverOverlapSize);
308  overlappingMatrix_ = new OverlappingMatrix(M,
309  borderListCreator.borderList(),
310  borderListCreator.blackList(),
311  overlapSize);
312 
313  // create the overlapping vectors for the residual and the
314  // solution
315  overlappingb_ = new OverlappingVector(overlappingMatrix_->overlap());
316  overlappingx_ = new OverlappingVector(*overlappingb_);
317 
318  // writeOverlapToVTK_();
319  }
320 
321  void rescale_()
322  {
323  const auto& overlap = overlappingMatrix_->overlap();
324  for (unsigned domesticRowIdx = 0; domesticRowIdx < overlap.numLocal(); ++domesticRowIdx) {
325  Index nativeRowIdx = overlap.domesticToNative(static_cast<Index>(domesticRowIdx));
326  auto& row = (*overlappingMatrix_)[domesticRowIdx];
327 
328  auto colIt = row.begin();
329  const auto& colEndIt = row.end();
330  for (; colIt != colEndIt; ++ colIt) {
331  auto& entry = *colIt;
332  for (unsigned i = 0; i < entry.rows; ++i)
333  entry[i] *= simulator_.model().eqWeight(nativeRowIdx, i);
334  }
335 
336  auto& rhsEntry = (*overlappingb_)[domesticRowIdx];
337  for (unsigned i = 0; i < rhsEntry.size(); ++i)
338  rhsEntry[i] *= simulator_.model().eqWeight(nativeRowIdx, i);
339  }
340  }
341 
342  void cleanup_()
343  {
344  // create the overlapping Jacobian matrix and vectors
345  delete overlappingMatrix_;
346  delete overlappingb_;
347  delete overlappingx_;
348 
349  overlappingMatrix_ = 0;
350  overlappingb_ = 0;
351  overlappingx_ = 0;
352  }
353 
354  std::shared_ptr<ParallelPreconditioner> preparePreconditioner_()
355  {
356  int preconditionerIsReady = 1;
357  try {
358  // update sequential preconditioner
359  precWrapper_.prepare(*overlappingMatrix_);
360  }
361  catch (const Dune::Exception& e) {
362  std::cout << "Preconditioner threw exception \"" << e.what()
363  << " on rank " << overlappingMatrix_->overlap().myRank()
364  << "\n" << std::flush;
365  preconditionerIsReady = 0;
366  }
367 
368  // make sure that the preconditioner is also ready on all peer
369  // ranks.
370  preconditionerIsReady = simulator_.gridView().comm().min(preconditionerIsReady);
371  if (!preconditionerIsReady)
372  OPM_THROW(Opm::NumericalProblem, "Creating the preconditioner failed");
373 
374  // create the parallel preconditioner
375  return std::make_shared<ParallelPreconditioner>(precWrapper_.get(), overlappingMatrix_->overlap());
376  }
377 
378  void cleanupPreconditioner_()
379  {
380  precWrapper_.cleanup();
381  }
382 
383  void writeOverlapToVTK_()
384  {
385  for (int lookedAtRank = 0;
386  lookedAtRank < simulator_.gridView().comm().size(); ++lookedAtRank) {
387  std::cout << "writing overlap for rank " << lookedAtRank << "\n" << std::flush;
388  typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > VtkField;
389  int n = simulator_.gridView().size(/*codim=*/dimWorld);
390  VtkField isInOverlap(n);
391  VtkField rankField(n);
392  isInOverlap = 0.0;
393  rankField = 0.0;
394  assert(rankField.two_norm() == 0.0);
395  assert(isInOverlap.two_norm() == 0.0);
396  auto vIt = simulator_.gridView().template begin</*codim=*/dimWorld>();
397  const auto& vEndIt = simulator_.gridView().template end</*codim=*/dimWorld>();
398  const auto& overlap = overlappingMatrix_->overlap();
399  for (; vIt != vEndIt; ++vIt) {
400  int nativeIdx = simulator_.model().vertexMapper().map(*vIt);
401  int localIdx = overlap.foreignOverlap().nativeToLocal(nativeIdx);
402  if (localIdx < 0)
403  continue;
404  rankField[nativeIdx] = simulator_.gridView().comm().rank();
405  if (overlap.peerHasIndex(lookedAtRank, localIdx))
406  isInOverlap[nativeIdx] = 1.0;
407  }
408 
409  typedef Dune::VTKWriter<GridView> VtkWriter;
410  VtkWriter writer(simulator_.gridView(), Dune::VTK::conforming);
411  writer.addVertexData(isInOverlap, "overlap");
412  writer.addVertexData(rankField, "rank");
413 
414  std::ostringstream oss;
415  oss << "overlap_rank=" << lookedAtRank;
416  writer.write(oss.str().c_str(), Dune::VTK::ascii);
417  }
418  }
419 
420  const Simulator& simulator_;
421  int gridSequenceNumber_;
422 
423  OverlappingMatrix *overlappingMatrix_;
424  OverlappingVector *overlappingb_;
425  OverlappingVector *overlappingx_;
426 
427  PreconditionerWrapper precWrapper_;
428 };
429 }} // namespace Linear, Ewoms
430 
431 namespace Ewoms {
432 namespace Properties {
434 SET_INT_PROP(ParallelBaseLinearSolver, LinearSolverVerbosity, 0);
435 
437 SET_SCALAR_PROP(ParallelBaseLinearSolver, PreconditionerRelaxation, 1.0);
438 
440 SET_INT_PROP(ParallelBaseLinearSolver, PreconditionerOrder, 0);
441 
444 SET_TYPE_PROP(ParallelBaseLinearSolver,
445  LinearSolverScalar,
446  typename GET_PROP_TYPE(TypeTag, Scalar));
447 
448 SET_PROP(ParallelBaseLinearSolver, OverlappingMatrix)
449 {
450  static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq);
451  typedef typename GET_PROP_TYPE(TypeTag, LinearSolverScalar) LinearSolverScalar;
452  typedef Dune::FieldMatrix<LinearSolverScalar, numEq, numEq> MatrixBlock;
453  typedef Dune::BCRSMatrix<MatrixBlock> NonOverlappingMatrix;
455 };
456 
457 SET_TYPE_PROP(ParallelBaseLinearSolver,
458  Overlap,
459  typename GET_PROP_TYPE(TypeTag, OverlappingMatrix)::Overlap);
460 
461 SET_PROP(ParallelBaseLinearSolver, OverlappingVector)
462 {
463  static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq);
464  typedef typename GET_PROP_TYPE(TypeTag, LinearSolverScalar) LinearSolverScalar;
465  typedef Dune::FieldVector<LinearSolverScalar, numEq> VectorBlock;
466  typedef typename GET_PROP_TYPE(TypeTag, Overlap) Overlap;
468 };
469 
470 SET_PROP(ParallelBaseLinearSolver, OverlappingScalarProduct)
471 {
472  typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
473  typedef typename GET_PROP_TYPE(TypeTag, Overlap) Overlap;
475 };
476 
477 SET_PROP(ParallelBaseLinearSolver, OverlappingLinearOperator)
478 {
479  typedef typename GET_PROP_TYPE(TypeTag, OverlappingMatrix) OverlappingMatrix;
480  typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
481  typedef Ewoms::Linear::OverlappingOperator<OverlappingMatrix, OverlappingVector,
482  OverlappingVector> type;
483 };
484 
485 SET_TYPE_PROP(ParallelBaseLinearSolver,
486  PreconditionerWrapper,
487  Ewoms::Linear::PreconditionerWrapperILU0<TypeTag>);
488 
490 SET_INT_PROP(ParallelBaseLinearSolver, LinearSolverOverlapSize, 2);
491 
493 SET_INT_PROP(ParallelBaseLinearSolver, LinearSolverMaxIterations, 1000);
494 } // namespace Properties
495 } // namespace Ewoms
496 
497 #endif
An overlap aware ISTL scalar product.
Definition: overlappingscalarproduct.hh:41
Definition: baseauxiliarymodule.hh:37
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
An overlap aware block vector.
Definition: overlappingblockvector.hh:49
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
A simple class which makes sure that a cleanup function is called once the object is destroyed...
An overlap aware linear operator usable by ISTL.
Definition: overlappingoperator.hh:39
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
An overlap aware block vector.
Provides the common code which is required by most linear solvers.
Definition: parallelbasebackend.hh:142
int Index
The type of an index of a degree of freedom.
Definition: overlaptypes.hh:44
#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
This file provides the infrastructure to retrieve run-time parameters.
bool solve(Vector &x)
Actually solve the linear system of equations.
Definition: parallelbasebackend.hh:248
An overlap aware block-compressed row storage (BCRS) matrix.
Definition: overlappingbcrsmatrix.hh:54
An overlap aware linear operator usable by ISTL.
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition: simulator.hh:183
void eraseMatrix()
Causes the solve() method to discared the structure of the linear system of equations the next time i...
Definition: parallelbasebackend.hh:205
Provides wrapper classes for the (non-AMG) preconditioners provided by dune-istl. ...
#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
An overlap aware preconditioner for any ISTL linear solver.
Definition: overlappingpreconditioner.hh:44
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Provides the magic behind the eWoms property system.
An overlap aware preconditioner for any ISTL 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
A simple class which makes sure that a cleanup function is called once the object is destroyed...
Definition: genericguard.hh:42
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
An overlap aware block-compressed row storage (BCRS) matrix.
#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.
An overlap aware ISTL scalar product.