26 #ifndef EWOMS_PARALLEL_ITERATIVE_BACKEND_HH
27 #define EWOMS_PARALLEL_ITERATIVE_BACKEND_HH
41 #include <dune/grid/io/file/vtk/vtkwriter.hh>
43 #include <dune/istl/preconditioners.hh>
45 #include <dune/common/shared_ptr.hh>
46 #include <dune/common/fvector.hh>
52 namespace Properties {
161 template <
class TypeTag>
164 typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) Implementation;
167 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
168 typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) Matrix;
169 typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
170 typedef typename GET_PROP_TYPE(TypeTag, BorderListCreator) BorderListCreator;
171 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
173 typedef typename GET_PROP_TYPE(TypeTag, Overlap) Overlap;
174 typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
175 typedef typename GET_PROP_TYPE(TypeTag, OverlappingMatrix) OverlappingMatrix;
177 typedef typename GET_PROP_TYPE(TypeTag, PreconditionerWrapper) PreconditionerWrapper;
178 typedef typename PreconditionerWrapper::SequentialPreconditioner SequentialPreconditioner;
180 typedef typename GET_PROP_TYPE(TypeTag, LinearSolverWrapper) LinearSolverWrapper;
190 enum { dimWorld = GridView::dimensionworld };
194 : simulator_(simulator)
196 overlappingMatrix_ = 0;
210 "The maximum allowed error between of the linear solver");
212 "The size of the algebraic overlap for the linear solver");
214 "The maximum number of iterations of the linear solver");
216 "The verbosity level of the linear solver");
218 LinearSolverWrapper::registerParameters();
219 PreconditionerWrapper::registerParameters();
241 bool solve(
const Matrix &M, Vector &x, Vector &b)
243 Scalar oldSingularLimit = Dune::FMatrixPrecision<Scalar>::singular_limit();
244 Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-50);
246 if (!overlappingMatrix_) {
255 overlappingMatrix_->assignAdd(M);
256 overlappingb_->assignAddBorder(b);
262 overlappingb_->assignTo(b);
264 (*overlappingx_) = 0.0;
266 int preconditionerIsReady = 1;
269 precWrapper_.prepare(*overlappingMatrix_);
271 catch (
const Dune::Exception &e) {
272 std::cout <<
"Preconditioner threw exception \"" << e.what()
273 <<
" on rank " << overlappingMatrix_->overlap().myRank()
274 <<
"\n" << std::flush;
275 preconditionerIsReady = 0;
280 preconditionerIsReady = simulator_.
gridView().comm().min(preconditionerIsReady);
281 if (!preconditionerIsReady) {
282 Dune::FMatrixPrecision<Scalar>::set_singular_limit(oldSingularLimit);
287 ParallelPreconditioner parPreCond(precWrapper_.get(),
288 overlappingMatrix_->overlap());
291 ParallelScalarProduct parScalarProduct(overlappingMatrix_->overlap());
292 ParallelOperator parOperator(*overlappingMatrix_);
295 auto &solver = solverWrapper_.get(parOperator, parScalarProduct, parPreCond);
301 OverlappingVector residWeightVec(*overlappingx_);
302 residWeightVec = 0.0;
303 const auto &overlap = overlappingMatrix_->overlap();
304 for (
unsigned localIdx = 0; localIdx < unsigned(overlap.numLocal()); ++localIdx) {
305 int nativeIdx = overlap.domesticToNative(localIdx);
306 for (
int eqIdx = 0; eqIdx < Vector::block_type::dimension; ++eqIdx) {
307 residWeightVec[localIdx][eqIdx] =
308 this->simulator_.
model().eqWeight(nativeIdx, eqIdx);
312 Scalar linearSolverTolerance =
EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance);
313 Scalar linearSolverAbsTolerance = simulator_.
model().newtonMethod().tolerance() / 100.0;
314 Scalar linearSolverFixPointTolerance = 100*std::numeric_limits<Scalar>::epsilon();
315 typedef typename GridView::CollectiveCommunication Comm;
320 linearSolverFixPointTolerance,
321 linearSolverTolerance,
322 linearSolverAbsTolerance);
329 solver.setConvergenceCriterion(std::shared_ptr<ConvergenceCriterion>(convCrit));
332 Dune::InverseOperatorResult result;
333 int solverSucceeded = 1;
335 solver.apply(*overlappingx_, *overlappingb_, result);
336 solverSucceeded = simulator_.
gridView().comm().min(solverSucceeded);
338 catch (
const Dune::Exception &) {
340 solverSucceeded = simulator_.
gridView().comm().min(solverSucceeded);
345 solverWrapper_.cleanup();
346 precWrapper_.cleanup();
348 if (!solverSucceeded) {
349 Dune::FMatrixPrecision<Scalar>::set_singular_limit(oldSingularLimit);
354 overlappingx_->assignTo(x);
358 Dune::FMatrixPrecision<Scalar>::set_singular_limit(oldSingularLimit);
361 return result.converged;
365 Implementation &asImp_()
366 {
return *
static_cast<Implementation *
>(
this); }
368 const Implementation &asImp_()
const
369 {
return *
static_cast<const Implementation *
>(
this); }
371 void prepare_(
const Matrix &M)
373 BorderListCreator borderListCreator(simulator_.
gridView(),
374 simulator_.
model().dofMapper());
377 int overlapSize =
EWOMS_GET_PARAM(TypeTag,
int, LinearSolverOverlapSize);
378 overlappingMatrix_ =
new OverlappingMatrix(M,
379 borderListCreator.borderList(),
380 borderListCreator.blackList(),
385 overlappingb_ =
new OverlappingVector(overlappingMatrix_->overlap());
386 overlappingx_ =
new OverlappingVector(*overlappingb_);
394 delete overlappingMatrix_;
395 delete overlappingb_;
396 delete overlappingx_;
398 overlappingMatrix_ = 0;
403 void writeOverlapToVTK_()
405 for (
int lookedAtRank = 0;
406 lookedAtRank < simulator_.
gridView().comm().size(); ++lookedAtRank) {
407 std::cout <<
"writing overlap for rank " << lookedAtRank <<
"\n" << std::flush;
408 typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > VtkField;
409 int n = simulator_.
gridView().size(dimWorld);
410 VtkField isInOverlap(n);
411 VtkField rankField(n);
414 assert(rankField.two_norm() == 0.0);
415 assert(isInOverlap.two_norm() == 0.0);
416 auto vIt = simulator_.
gridView().template begin<dimWorld>();
417 const auto &vEndIt = simulator_.
gridView().template end<dimWorld>();
418 const auto &overlap = overlappingMatrix_->overlap();
419 for (; vIt != vEndIt; ++vIt) {
420 int nativeIdx = simulator_.
model().vertexMapper().map(*vIt);
421 int localIdx = overlap.foreignOverlap().nativeToLocal(nativeIdx);
424 rankField[nativeIdx] = simulator_.
gridView().comm().rank();
425 if (overlap.peerHasIndex(lookedAtRank, localIdx))
426 isInOverlap[nativeIdx] = 1.0;
429 typedef Dune::VTKWriter<GridView> VtkWriter;
430 VtkWriter writer(simulator_.
gridView(), Dune::VTK::conforming);
431 writer.addVertexData(isInOverlap,
"overlap");
432 writer.addVertexData(rankField,
"rank");
434 std::ostringstream oss;
435 oss <<
"overlap_rank=" << lookedAtRank;
436 writer.write(oss.str().c_str(), Dune::VTK::ascii);
440 const Simulator &simulator_;
442 OverlappingMatrix *overlappingMatrix_;
443 OverlappingVector *overlappingb_;
444 OverlappingVector *overlappingx_;
446 PreconditionerWrapper precWrapper_;
447 LinearSolverWrapper solverWrapper_;
453 #define EWOMS_WRAP_ISTL_SOLVER(SOLVER_NAME, ISTL_SOLVER_NAME) \
454 template <class TypeTag> \
455 class SolverWrapper##SOLVER_NAME \
457 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; \
458 typedef typename GET_PROP_TYPE(TypeTag, \
459 OverlappingMatrix) OverlappingMatrix; \
460 typedef typename GET_PROP_TYPE(TypeTag, \
461 OverlappingVector) OverlappingVector; \
463 typedef ISTL_SOLVER_NAME<OverlappingVector> ParallelSolver; \
466 SolverWrapper##SOLVER_NAME() \
469 static void registerParameters() \
472 template <class LinearOperator, class ScalarProduct, class Preconditioner> \
473 ParallelSolver &get(LinearOperator &parOperator, \
474 ScalarProduct &parScalarProduct, \
475 Preconditioner &parPreCond) \
477 Scalar tolerance = EWOMS_GET_PARAM(TypeTag, Scalar, \
478 LinearSolverTolerance); \
479 int maxIter = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIterations);\
482 if (parOperator.overlap().myRank() == 0) \
483 verbosity = EWOMS_GET_PARAM(TypeTag, int, LinearSolverVerbosity); \
484 solver_ = new ParallelSolver(parOperator, parScalarProduct, \
485 parPreCond, tolerance, maxIter, \
492 { delete solver_; } \
495 ParallelSolver *solver_; \
505 #undef EWOMS_WRAP_ISTL_SOLVER
506 #undef EWOMS_ISTL_SOLVER_TYPDEF
508 #define EWOMS_WRAP_ISTL_PRECONDITIONER(PREC_NAME, ISTL_PREC_TYPE) \
509 template <class TypeTag> \
510 class PreconditionerWrapper##PREC_NAME \
512 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; \
513 typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix; \
514 typedef typename GET_PROP_TYPE(TypeTag, \
515 OverlappingVector) OverlappingVector; \
518 typedef ISTL_PREC_TYPE<JacobianMatrix, OverlappingVector, \
519 OverlappingVector> SequentialPreconditioner; \
520 PreconditionerWrapper##PREC_NAME() \
523 static void registerParameters() \
525 EWOMS_REGISTER_PARAM(TypeTag, int, PreconditionerOrder, \
526 "The order of the preconditioner"); \
527 EWOMS_REGISTER_PARAM(TypeTag, Scalar, PreconditionerRelaxation, \
528 "The relaxation factor of the " \
532 void prepare(JacobianMatrix &matrix) \
534 int order = EWOMS_GET_PARAM(TypeTag, int, PreconditionerOrder); \
535 Scalar relaxationFactor = EWOMS_GET_PARAM(TypeTag, Scalar, PreconditionerRelaxation); \
536 seqPreCond_ = new SequentialPreconditioner(matrix, order, \
540 SequentialPreconditioner &get() \
541 { return *seqPreCond_; } \
544 { delete seqPreCond_; } \
547 SequentialPreconditioner *seqPreCond_; \
552 #define EWOMS_WRAP_ISTL_SIMPLE_PRECONDITIONER(PREC_NAME, ISTL_PREC_TYPE) \
553 template <class TypeTag> \
554 class PreconditionerWrapper##PREC_NAME \
556 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; \
557 typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix; \
558 typedef typename GET_PROP_TYPE(TypeTag, \
559 OverlappingVector) OverlappingVector; \
562 typedef ISTL_PREC_TYPE<JacobianMatrix, OverlappingVector, \
563 OverlappingVector> SequentialPreconditioner; \
564 PreconditionerWrapper##PREC_NAME() \
567 static void registerParameters() \
569 EWOMS_REGISTER_PARAM(TypeTag, Scalar, PreconditionerRelaxation, \
570 "The relaxation factor of the " \
574 void prepare(JacobianMatrix &matrix) \
576 Scalar relaxationFactor= EWOMS_GET_PARAM(TypeTag, Scalar, PreconditionerRelaxation);\
577 seqPreCond_ = new SequentialPreconditioner(matrix, \
581 SequentialPreconditioner &get() \
582 { return *seqPreCond_; } \
585 { delete seqPreCond_; } \
588 SequentialPreconditioner *seqPreCond_; \
600 #undef EWOMS_WRAP_ISTL_PRECONDITIONER
605 namespace Properties {
607 SET_INT_PROP(ParallelIterativeLinearSolver, LinearSolverVerbosity, 0);
610 SET_SCALAR_PROP(ParallelIterativeLinearSolver, PreconditionerRelaxation, 1.0);
613 SET_INT_PROP(ParallelIterativeLinearSolver, PreconditionerOrder, 0);
616 SET_INT_PROP(ParallelIterativeLinearSolver, GMResRestart, 10);
618 SET_TYPE_PROP(ParallelIterativeLinearSolver, OverlappingMatrix,
620 TypeTag, JacobianMatrix)>);
622 typename GET_PROP_TYPE(TypeTag, OverlappingMatrix)::Overlap);
623 SET_PROP(ParallelIterativeLinearSolver, OverlappingVector)
625 typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
630 SET_PROP(ParallelIterativeLinearSolver, OverlappingScalarProduct)
632 typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
636 SET_PROP(ParallelIterativeLinearSolver, OverlappingLinearOperator)
638 typedef typename GET_PROP_TYPE(TypeTag, OverlappingMatrix) OverlappingMatrix;
639 typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
641 OverlappingVector> type;
644 SET_TYPE_PROP(ParallelIterativeLinearSolver, LinearSolverBackend,
646 SET_TYPE_PROP(ParallelIterativeLinearSolver, LinearSolverWrapper,
647 Ewoms::Linear::SolverWrapperBiCGStab<TypeTag>);
648 SET_TYPE_PROP(ParallelIterativeLinearSolver, PreconditionerWrapper,
649 Ewoms::Linear::PreconditionerWrapperILU0<TypeTag>);
652 SET_SCALAR_PROP(ParallelIterativeLinearSolver, LinearSolverOverlapSize, 2);
655 SET_INT_PROP(ParallelIterativeLinearSolver, LinearSolverMaxIterations, 250);
#define EWOMS_WRAP_ISTL_SIMPLE_PRECONDITIONER(PREC_NAME, ISTL_PREC_TYPE)
Definition: paralleliterativebackend.hh:552
An overlap aware ISTL scalar product.
Definition: overlappingscalarproduct.hh:42
An ISTL preconditioner that solves the linear system of equations locally on each rank...
An overlap aware block-compressed row storage (BCRS) matrix.
An overlap aware linear operator usable by ISTL.
SET_SCALAR_PROP(NumericModel, EndTime,-1e100)
The default value for the simulation's end time.
Base class for all convergence criteria which only defines an virtual API.
Definition: convergencecriterion.hh:51
An overlap aware preconditioner for any ISTL linear solver.
Definition: overlappingpreconditioner.hh:42
Copy of dune-istl's linear solvers with added support for pluggable convergence criteria.
An overlap aware block-compressed row storage (BCRS) matrix.
Definition: overlappingbcrsmatrix.hh:52
void setStructureMatrix(const Matrix &M)
Set the structure of the linear system of equations to be solved.
Definition: paralleliterativebackend.hh:230
An overlap aware block vector.
Definition: overlappingblockvector.hh:47
#define EWOMS_WRAP_ISTL_SOLVER(SOLVER_NAME, ISTL_SOLVER_NAME)
Macro to create a wrapper around an ISTL solver.
Definition: paralleliterativebackend.hh:453
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:485
ParallelIterativeSolverBackend(const Simulator &simulator)
Definition: paralleliterativebackend.hh:193
SET_PROP(NumericModel, ParameterTree)
Set the ParameterTree property.
Definition: basicproperties.hh:117
An overlap aware ISTL scalar product.
SET_INT_PROP(NumericModel, GridGlobalRefinements, 0)
Convergence criterion which looks at the weighted absolute value of the residual. ...
Definition: weightedresidreductioncriterion.hh:54
Definition: cartesianindexmapper.hh:31
NEW_PROP_TAG(Grid)
The type of the DUNE grid.
bool solve(const Matrix &M, Vector &x, Vector &b)
Actually solve the linear system of equations.
Definition: paralleliterativebackend.hh:241
This file provides the infrastructure to retrieve run-time parameters.
SET_TYPE_PROP(NumericModel, Scalar, double)
Set the default type of scalar values to double.
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:73
Definition: baseauxiliarymodule.hh:35
Preconditioned loop solver.
Definition: solvers.hh:167
Implements a generic linear solver abstraction.
Definition: paralleliterativebackend.hh:162
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:64
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:176
~ParallelIterativeSolverBackend()
Definition: paralleliterativebackend.hh:201
static void registerParameters()
Register all run-time parameters for the linear solver.
Definition: paralleliterativebackend.hh:207
GridView & gridView()
Return the grid view for which the simulation is done.
Definition: simulator.hh:164
Provides the magic behind the eWoms property system.
An overlap aware block vector.
#define EWOMS_WRAP_ISTL_PRECONDITIONER(PREC_NAME, ISTL_PREC_TYPE)
Definition: paralleliterativebackend.hh:508
An overlap aware preconditioner for any ISTL linear solver.
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:95
An overlap aware linear operator usable by ISTL.
Definition: overlappingoperator.hh:37