27#ifndef EWOMS_PARALLEL_BASE_BACKEND_HH
28#define EWOMS_PARALLEL_BASE_BACKEND_HH
30#include <dune/common/fvector.hh>
31#include <dune/common/version.hh>
33#include <dune/grid/io/file/vtk/vtkwriter.hh>
35#include <opm/common/Exceptions.hpp>
45#include <opm/simulators/linalg/matrixblock.hh>
56namespace Opm::Properties {
64template<
class TypeTag>
107template <
class TypeTag>
126 using SequentialPreconditioner =
typename PreconditionerWrapper::SequentialPreconditioner;
134 enum { dimWorld = GridView::dimensionworld };
138 : simulator_(simulator)
139 , gridSequenceNumber_( -1 )
140 , lastIterations_( -1 )
142 overlappingMatrix_ =
nullptr;
143 overlappingb_ =
nullptr;
144 overlappingx_ =
nullptr;
155 Parameters::Register<Parameters::LinearSolverTolerance<Scalar>>
156 (
"The maximum allowed error between of the linear solver");
157 Parameters::Register<Parameters::LinearSolverAbsTolerance<Scalar>>
158 (
"The maximum accepted error of the norm of the residual");
159 Parameters::Register<Parameters::LinearSolverOverlapSize>
160 (
"The size of the algebraic overlap for the linear solver");
161 Parameters::Register<Parameters::LinearSolverMaxIterations>
162 (
"The maximum number of iterations of the linear solver");
163 Parameters::Register<Parameters::LinearSolverVerbosity>
164 (
"The verbosity level of the linear solver");
166 PreconditionerWrapper::registerParameters();
182 void prepare(
const SparseMatrixAdapter& M,
const Vector& )
185 int curSeqNum = simulator_.vanguard().gridSequenceNumber();
186 if (gridSequenceNumber_ == curSeqNum && overlappingMatrix_)
192 gridSequenceNumber_ = curSeqNum;
194 BorderListCreator borderListCreator(simulator_.gridView(),
195 simulator_.model().dofMapper());
198 unsigned overlapSize = Parameters::Get<Parameters::LinearSolverOverlapSize>();
199 overlappingMatrix_ =
new OverlappingMatrix(M.istlMatrix(),
200 borderListCreator.borderList(),
201 borderListCreator.blackList(),
206 overlappingb_ =
new OverlappingVector(overlappingMatrix_->overlap());
207 overlappingx_ =
new OverlappingVector(*overlappingb_);
221 overlappingb_->assignAddBorder(b);
232 overlappingb_->assignTo(b);
243 overlappingMatrix_->assignFromNative(M.istlMatrix());
244 overlappingMatrix_->syncAdd();
254 (*overlappingx_) = 0.0;
256 auto parPreCond = asImp_().preparePreconditioner_();
257 auto precondCleanupFn = [
this]() ->
void
258 { this->asImp_().cleanupPreconditioner_(); };
259 auto precondCleanupGuard = Opm::make_guard(precondCleanupFn);
265 auto solver = asImp_().prepareSolver_(parOperator,
269 auto cleanupSolverFn =
271 { this->asImp_().cleanupSolver_(); };
272 GenericGuard<
decltype(cleanupSolverFn)> solverGuard(cleanupSolverFn);
275 auto result = asImp_().runSolver_(solver);
277 lastIterations_ = result.second;
280 overlappingx_->assignTo(x);
290 {
return lastIterations_; }
293 Implementation& asImp_()
294 {
return *
static_cast<Implementation *
>(
this); }
296 const Implementation& asImp_()
const
297 {
return *
static_cast<const Implementation *
>(
this); }
302 delete overlappingMatrix_;
303 delete overlappingb_;
304 delete overlappingx_;
306 overlappingMatrix_ = 0;
311 std::shared_ptr<ParallelPreconditioner> preparePreconditioner_()
313 int preconditionerIsReady = 1;
316 precWrapper_.prepare(*overlappingMatrix_);
318 catch (
const Dune::Exception& e) {
319 std::cout <<
"Preconditioner threw exception \"" << e.what()
320 <<
" on rank " << overlappingMatrix_->overlap().myRank()
321 <<
"\n" << std::flush;
322 preconditionerIsReady = 0;
327 preconditionerIsReady = simulator_.gridView().comm().min(preconditionerIsReady);
328 if (!preconditionerIsReady)
329 throw NumericalProblem(
"Creating the preconditioner failed");
332 return std::make_shared<ParallelPreconditioner>(precWrapper_.get(), overlappingMatrix_->overlap());
335 void cleanupPreconditioner_()
337 precWrapper_.cleanup();
340 void writeOverlapToVTK_()
342 for (
int lookedAtRank = 0;
343 lookedAtRank < simulator_.gridView().comm().size(); ++lookedAtRank) {
344 std::cout <<
"writing overlap for rank " << lookedAtRank <<
"\n" << std::flush;
345 using VtkField = Dune::BlockVector<Dune::FieldVector<Scalar, 1> >;
346 int n = simulator_.gridView().size(dimWorld);
347 VtkField isInOverlap(n);
348 VtkField rankField(n);
351 assert(rankField.two_norm() == 0.0);
352 assert(isInOverlap.two_norm() == 0.0);
353 auto vIt = simulator_.gridView().template begin<dimWorld>();
354 const auto& vEndIt = simulator_.gridView().template end<dimWorld>();
355 const auto& overlap = overlappingMatrix_->overlap();
356 for (; vIt != vEndIt; ++vIt) {
357 int nativeIdx = simulator_.model().vertexMapper().map(*vIt);
358 int localIdx = overlap.foreignOverlap().nativeToLocal(nativeIdx);
361 rankField[nativeIdx] = simulator_.gridView().comm().rank();
362 if (overlap.peerHasIndex(lookedAtRank, localIdx))
363 isInOverlap[nativeIdx] = 1.0;
366 using VtkWriter = Dune::VTKWriter<GridView>;
367 VtkWriter writer(simulator_.gridView(), Dune::VTK::conforming);
368 writer.addVertexData(isInOverlap,
"overlap");
369 writer.addVertexData(rankField,
"rank");
371 std::ostringstream oss;
372 oss <<
"overlap_rank=" << lookedAtRank;
373 writer.write(oss.str().c_str(), Dune::VTK::ascii);
377 const Simulator& simulator_;
378 int gridSequenceNumber_;
379 size_t lastIterations_;
381 OverlappingMatrix *overlappingMatrix_;
382 OverlappingVector *overlappingb_;
383 OverlappingVector *overlappingx_;
385 PreconditionerWrapper precWrapper_;
389namespace Opm::Properties {
393template<
class TypeTag>
397template<
class TypeTag>
404 using NonOverlappingMatrix = Dune::BCRSMatrix<MatrixBlock>;
410template<
class TypeTag>
411struct Overlap<TypeTag, TTag::ParallelBaseLinearSolver>
414template<
class TypeTag>
419 using VectorBlock = Dune::FieldVector<LinearSolverScalar, numEq>;
424template<
class TypeTag>
432template<
class TypeTag>
441template<
class TypeTag>
A simple class which makes sure that a cleanup function is called once the object is destroyed.
Definition genericguard.hh:43
A sparse matrix interface backend for BCRSMatrix from dune-istl.
Definition istlsparsematrixadapter.hh:43
An overlap aware block-compressed row storage (BCRS) matrix.
Definition overlappingbcrsmatrix.hh:54
An overlap aware block vector.
Definition overlappingblockvector.hh:50
An overlap aware linear operator usable by ISTL.
Definition overlappingoperator.hh:42
An overlap aware preconditioner for any ISTL linear solver.
Definition overlappingpreconditioner.hh:48
An overlap aware ISTL scalar product.
Definition overlappingscalarproduct.hh:42
Provides the common code which is required by most linear solvers.
Definition parallelbasebackend.hh:109
void getResidual(Vector &b) const
Retrieve the synchronized internal residual vector.
Definition parallelbasebackend.hh:229
size_t iterations() const
Return number of iterations used during last solve.
Definition parallelbasebackend.hh:289
void eraseMatrix()
Causes the solve() method to discared the structure of the linear system of equations the next time i...
Definition parallelbasebackend.hh:173
void setResidual(const Vector &b)
Assign values to the internal data structure for the residual vector.
Definition parallelbasebackend.hh:217
void prepare(const SparseMatrixAdapter &M, const Vector &)
Set up the internal data structures required for the linear solver.
Definition parallelbasebackend.hh:182
void setMatrix(const SparseMatrixAdapter &M)
Sets the values of the residual's Jacobian matrix.
Definition parallelbasebackend.hh:241
static void registerParameters()
Register all run-time parameters for the linear solver.
Definition parallelbasebackend.hh:153
bool solve(Vector &x)
Actually solve the linear system of equations.
Definition parallelbasebackend.hh:252
Definition istlpreconditionerwrappers.hh:153
Definition MatrixMarketSpecializations.hpp:17
A simple class which makes sure that a cleanup function is called once the object is destroyed.
Provides wrapper classes for the (non-AMG) preconditioners provided by dune-istl.
A sparse matrix interface backend for BCRSMatrix from dune-istl.
Declares the parameters for the black oil model.
Declares the properties required by the black oil model.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
An overlap aware block-compressed row storage (BCRS) matrix.
An overlap aware block vector.
An overlap aware linear operator usable by ISTL.
An overlap aware preconditioner for any ISTL linear solver.
An overlap aware ISTL scalar product.
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
The floating point type used internally by the linear solver.
Definition linalgproperties.hh:46
Definition linalgproperties.hh:60
Definition linalgproperties.hh:63
Definition linalgproperties.hh:66
Definition linalgproperties.hh:69
Definition linalgproperties.hh:72
the preconditioner used by the linear solver
Definition linalgproperties.hh:42
The class that allows to manipulate sparse matrices.
Definition linalgproperties.hh:50
Definition parallelbasebackend.hh:59