27 #ifndef EWOMS_PARALLEL_BASE_BACKEND_HH 28 #define EWOMS_PARALLEL_BASE_BACKEND_HH 42 #include <dune/grid/io/file/vtk/vtkwriter.hh> 44 #include <dune/common/fvector.hh> 51 namespace Properties {
141 template <
class TypeTag>
145 typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) Implementation;
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;
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;
158 typedef typename GET_PROP_TYPE(TypeTag, PreconditionerWrapper) PreconditionerWrapper;
159 typedef typename PreconditionerWrapper::SequentialPreconditioner SequentialPreconditioner;
169 enum { dimWorld = GridView::dimensionworld };
173 : simulator_(simulator)
174 , gridSequenceNumber_( -1 )
176 overlappingMatrix_ =
nullptr;
177 overlappingb_ =
nullptr;
178 overlappingx_ =
nullptr;
190 "The maximum allowed error between of the linear solver");
192 "The size of the algebraic overlap for the linear solver");
194 "The maximum number of iterations of the linear solver");
196 "The verbosity level of the linear solver");
198 PreconditionerWrapper::registerParameters();
208 void prepareMatrix(
const Matrix& M)
217 overlappingMatrix_->assignFromNative(M);
223 overlappingMatrix_->syncAdd();
225 overlappingb_->sync();
228 void prepareRhs(
const Matrix& M, Vector& b)
234 overlappingb_->assignAddBorder(b);
240 overlappingb_->assignTo(b);
250 (*overlappingx_) = 0.0;
252 auto parPreCond = asImp_().preparePreconditioner_();
254 auto cleanupPrecondFn =
256 { this->asImp_().cleanupPreconditioner_(); };
265 auto solver = asImp_().prepareSolver_(parOperator,
269 auto cleanupSolverFn =
271 { this->asImp_().cleanupSolver_(); };
275 bool result = asImp_().runSolver_(solver);
278 overlappingx_->assignTo(x);
285 Implementation& asImp_()
286 {
return *
static_cast<Implementation *
>(
this); }
288 const Implementation& asImp_()
const 289 {
return *
static_cast<const Implementation *
>(
this); }
291 void prepare_(
const Matrix& M)
294 int curSeqNum = simulator_.
gridManager().gridSequenceNumber();
295 if( gridSequenceNumber_ == curSeqNum && overlappingMatrix_)
301 gridSequenceNumber_ = curSeqNum;
303 BorderListCreator borderListCreator(simulator_.
gridView(),
304 simulator_.
model().dofMapper());
307 unsigned overlapSize =
EWOMS_GET_PARAM(TypeTag,
unsigned, LinearSolverOverlapSize);
308 overlappingMatrix_ =
new OverlappingMatrix(M,
309 borderListCreator.borderList(),
310 borderListCreator.blackList(),
315 overlappingb_ =
new OverlappingVector(overlappingMatrix_->overlap());
316 overlappingx_ =
new OverlappingVector(*overlappingb_);
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];
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);
336 auto& rhsEntry = (*overlappingb_)[domesticRowIdx];
337 for (
unsigned i = 0; i < rhsEntry.size(); ++i)
338 rhsEntry[i] *= simulator_.
model().eqWeight(nativeRowIdx, i);
345 delete overlappingMatrix_;
346 delete overlappingb_;
347 delete overlappingx_;
349 overlappingMatrix_ = 0;
354 std::shared_ptr<ParallelPreconditioner> preparePreconditioner_()
356 int preconditionerIsReady = 1;
359 precWrapper_.prepare(*overlappingMatrix_);
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;
370 preconditionerIsReady = simulator_.
gridView().comm().min(preconditionerIsReady);
371 if (!preconditionerIsReady)
372 OPM_THROW(Opm::NumericalProblem,
"Creating the preconditioner failed");
375 return std::make_shared<ParallelPreconditioner>(precWrapper_.get(), overlappingMatrix_->overlap());
378 void cleanupPreconditioner_()
380 precWrapper_.cleanup();
383 void writeOverlapToVTK_()
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(dimWorld);
390 VtkField isInOverlap(n);
391 VtkField rankField(n);
394 assert(rankField.two_norm() == 0.0);
395 assert(isInOverlap.two_norm() == 0.0);
396 auto vIt = simulator_.
gridView().template begin<dimWorld>();
397 const auto& vEndIt = simulator_.
gridView().template end<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);
404 rankField[nativeIdx] = simulator_.
gridView().comm().rank();
405 if (overlap.peerHasIndex(lookedAtRank, localIdx))
406 isInOverlap[nativeIdx] = 1.0;
409 typedef Dune::VTKWriter<GridView> VtkWriter;
410 VtkWriter writer(simulator_.
gridView(), Dune::VTK::conforming);
411 writer.addVertexData(isInOverlap,
"overlap");
412 writer.addVertexData(rankField,
"rank");
414 std::ostringstream oss;
415 oss <<
"overlap_rank=" << lookedAtRank;
416 writer.write(oss.str().c_str(), Dune::VTK::ascii);
420 const Simulator& simulator_;
421 int gridSequenceNumber_;
423 OverlappingMatrix *overlappingMatrix_;
424 OverlappingVector *overlappingb_;
425 OverlappingVector *overlappingx_;
427 PreconditionerWrapper precWrapper_;
432 namespace Properties {
434 SET_INT_PROP(ParallelBaseLinearSolver, LinearSolverVerbosity, 0);
437 SET_SCALAR_PROP(ParallelBaseLinearSolver, PreconditionerRelaxation, 1.0);
440 SET_INT_PROP(ParallelBaseLinearSolver, PreconditionerOrder, 0);
448 SET_PROP(ParallelBaseLinearSolver, OverlappingMatrix)
451 typedef typename GET_PROP_TYPE(TypeTag, LinearSolverScalar) LinearSolverScalar;
452 typedef Dune::FieldMatrix<LinearSolverScalar, numEq, numEq> MatrixBlock;
453 typedef Dune::BCRSMatrix<MatrixBlock> NonOverlappingMatrix;
459 typename GET_PROP_TYPE(TypeTag, OverlappingMatrix)::Overlap);
461 SET_PROP(ParallelBaseLinearSolver, OverlappingVector)
464 typedef typename GET_PROP_TYPE(TypeTag, LinearSolverScalar) LinearSolverScalar;
465 typedef Dune::FieldVector<LinearSolverScalar, numEq> VectorBlock;
470 SET_PROP(ParallelBaseLinearSolver, OverlappingScalarProduct)
472 typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
477 SET_PROP(ParallelBaseLinearSolver, OverlappingLinearOperator)
479 typedef typename GET_PROP_TYPE(TypeTag, OverlappingMatrix) OverlappingMatrix;
480 typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
482 OverlappingVector> type;
486 PreconditionerWrapper,
487 Ewoms::Linear::PreconditionerWrapperILU0<TypeTag>);
490 SET_INT_PROP(ParallelBaseLinearSolver, LinearSolverOverlapSize, 2);
493 SET_INT_PROP(ParallelBaseLinearSolver, LinearSolverMaxIterations, 1000);
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.