25 #ifndef EWOMS_PARALLEL_AMG_BACKEND_HH
26 #define EWOMS_PARALLEL_AMG_BACKEND_HH
40 #include <dune/istl/preconditioners.hh>
41 #include <dune/istl/solvers.hh>
42 #include <dune/istl/paamg/amg.hh>
43 #include <dune/istl/paamg/pinfo.hh>
44 #include <dune/istl/schwarz.hh>
45 #include <dune/istl/owneroverlapcopy.hh>
47 #include <dune/common/parallel/indexset.hh>
48 #include <dune/common/version.hh>
49 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
50 #include <dune/common/parallel/mpicollectivecommunication.hh>
52 #include <dune/common/mpicollectivecommunication.hh>
59 template <
class TypeTag>
63 namespace Properties {
70 SET_INT_PROP(ParallelAmgLinearSolver, AmgCoarsenTarget, 5000);
83 template <
class TypeTag>
84 class ParallelAmgBackend
86 typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) Implementation;
89 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
90 typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) Matrix;
91 typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
92 typedef typename GET_PROP_TYPE(TypeTag, BorderListCreator) BorderListCreator;
93 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
95 typedef typename GET_PROP_TYPE(TypeTag, Overlap) Overlap;
96 typedef typename GET_PROP_TYPE(TypeTag, OverlappingVector) OverlappingVector;
97 typedef typename GET_PROP_TYPE(TypeTag, OverlappingMatrix) OverlappingMatrix;
102 typedef
Dune::SeqSOR<Matrix, Vector, Vector> SequentialSmoother;
109 typedef Dune::OwnerOverlapCopyCommunication<Ewoms::Linear::Index>
110 OwnerOverlapCopyCommunication;
111 typedef Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector,
112 OwnerOverlapCopyCommunication>
114 typedef Dune::OverlappingSchwarzScalarProduct<Vector,
115 OwnerOverlapCopyCommunication> FineScalarProduct;
116 typedef Dune::BlockPreconditioner<Vector,
118 OwnerOverlapCopyCommunication,
119 SequentialSmoother> ParallelSmoother;
120 typedef Dune::Amg::AMG<FineOperator,
123 OwnerOverlapCopyCommunication> AMG;
125 typedef Dune::MatrixAdapter<Matrix, Vector, Vector> FineOperator;
126 typedef Dune::SeqScalarProduct<Vector> FineScalarProduct;
127 typedef SequentialSmoother ParallelSmoother;
128 typedef Dune::Amg::AMG<FineOperator, Vector, ParallelSmoother> AMG;
133 : simulator_(simulator)
135 overlappingMatrix_ =
nullptr;
136 overlappingb_ =
nullptr;
137 overlappingx_ =
nullptr;
150 "The coarsening target for the agglomerations of "
151 "the AMG preconditioner");
173 bool solve(
const Matrix &M, Vector &x, Vector &b)
176 if (simulator_.
gridManager().gridView().comm().rank() == 0)
183 if (!overlappingMatrix_) {
192 overlappingMatrix_->assignAdd(M);
193 overlappingb_->assignAddBorder(b);
199 overlappingb_->assignTo(b);
201 (*overlappingx_) = 0.0;
208 FineOperator fineOperator(*overlappingMatrix_, *istlComm_);
210 FineOperator fineOperator(*overlappingMatrix_);
214 FineScalarProduct scalarProduct(*istlComm_);
216 FineScalarProduct scalarProduct;
219 setupAmg_(fineOperator);
224 if (verbosity > 1 && simulator_.
gridManager().gridView().comm().rank() == 0)
225 std::cout <<
"Creating the solver\n" << std::flush;
228 Scalar linearSolverTolerance =
EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance);
229 int maxIterations =
EWOMS_GET_PARAM(TypeTag,
int, LinearSolverMaxIterations);
230 SolverType solver(fineOperator,
233 linearSolverTolerance,
241 OverlappingVector residWeightVec(*overlappingx_);
242 residWeightVec = 0.0;
243 const auto &overlap = overlappingMatrix_->overlap();
244 for (
unsigned localIdx = 0; localIdx < unsigned(overlap.numLocal()); ++localIdx) {
245 int nativeIdx = overlap.domesticToNative(localIdx);
246 for (
int eqIdx = 0; eqIdx < Vector::block_type::dimension; ++eqIdx) {
247 residWeightVec[localIdx][eqIdx] =
248 this->simulator_.
model().eqWeight(nativeIdx, eqIdx);
252 Scalar linearSolverAbsTolerance = simulator_.
model().newtonMethod().tolerance() / 100.0;
253 Scalar linearSolverFixPointTolerance = 100*std::numeric_limits<Scalar>::epsilon();
254 typedef typename GridView::CollectiveCommunication Comm;
259 linearSolverFixPointTolerance,
260 linearSolverTolerance,
261 linearSolverAbsTolerance);
268 solver.setConvergenceCriterion(std::shared_ptr<ConvergenceCriterion>(convCrit));
270 Dune::InverseOperatorResult result;
271 int solverSucceeded = 1;
274 solver.apply(*overlappingx_, *overlappingb_, result);
275 solverSucceeded = simulator_.
gridManager().gridView().comm().min(solverSucceeded);
277 catch (
const Dune::Exception &)
280 solverSucceeded = simulator_.
gridManager().gridView().comm().min(solverSucceeded);
283 if (!solverSucceeded)
287 overlappingx_->assignTo(x);
289 return result.converged;
293 Implementation &asImp_()
294 {
return *
static_cast<Implementation *
>(
this); }
296 const Implementation &asImp_()
const
297 {
return *
static_cast<const Implementation *
>(
this); }
299 void prepare_(
const Matrix &M)
301 BorderListCreator borderListCreator(simulator_.
gridView(),
302 simulator_.
model().dofMapper());
304 auto &blackList = borderListCreator.blackList();
307 int overlapSize =
EWOMS_GET_PARAM(TypeTag,
int, LinearSolverOverlapSize);
308 overlappingMatrix_ =
new OverlappingMatrix(M,
309 borderListCreator.borderList(),
315 overlappingb_ =
new OverlappingVector(overlappingMatrix_->overlap());
316 overlappingx_ =
new OverlappingVector(*overlappingb_);
323 istlComm_ =
new OwnerOverlapCopyCommunication(MPI_COMM_WORLD);
324 setupAmgIndexSet_(overlappingMatrix_->overlap(), istlComm_->indexSet());
325 istlComm_->remoteIndices().template rebuild<false>();
332 delete overlappingMatrix_;
333 delete overlappingb_;
334 delete overlappingx_;
337 overlappingMatrix_ =
nullptr;
338 overlappingb_ =
nullptr;
339 overlappingx_ =
nullptr;
344 template <
class ParallelIndexSet>
345 void setupAmgIndexSet_(
const Overlap &overlap, ParallelIndexSet &istlIndices)
347 typedef Dune::OwnerOverlapCopyAttributeSet GridAttributes;
348 typedef Dune::OwnerOverlapCopyAttributeSet::AttributeSet GridAttributeSet;
351 istlIndices.beginResize();
352 for (
int curIdx = 0; curIdx < overlap.numDomestic(); ++curIdx) {
353 GridAttributeSet gridFlag = overlap.iAmMasterOf(curIdx)
354 ? GridAttributes::owner
355 : GridAttributes::copy;
359 bool isShared = overlap.isInOverlap(curIdx);
361 assert(curIdx ==
int(overlap.globalToDomestic(
362 overlap.domesticToGlobal(curIdx))));
363 istlIndices.add(overlap.domesticToGlobal(curIdx),
364 Dune::ParallelLocalIndex<GridAttributeSet>(
365 curIdx, gridFlag, isShared));
367 istlIndices.endResize();
371 void setupAmg_(FineOperator &fineOperator)
377 if (simulator_.
gridManager().gridView().comm().rank() == 0)
380 int rank = simulator_.
gridManager().gridView().comm().rank();
381 if (verbosity > 1 && rank == 0)
382 std::cout <<
"Setting up the AMG preconditioner\n";
384 typedef typename Dune::Amg::SmootherTraits<ParallelSmoother>::Arguments
387 SmootherArgs smootherArgs;
388 smootherArgs.iterations = 1;
389 smootherArgs.relaxationFactor = 1.0;
397 CoarsenCriterion<Dune::Amg::SymmetricCriterion<Matrix, Dune::Amg::FrobeniusNorm> >
400 CoarsenCriterion coarsenCriterion(15, coarsenTarget);
401 coarsenCriterion.setDefaultValuesAnisotropic(GridView::dimension,
404 coarsenCriterion.setDebugLevel(1);
406 coarsenCriterion.setDebugLevel(0);
409 coarsenCriterion.setMinCoarsenRate(1.05);
411 coarsenCriterion.setAccumulate(Dune::Amg::atOnceAccu);
412 coarsenCriterion.setSkipIsolated(
false);
416 amg_ =
new AMG(fineOperator, coarsenCriterion, smootherArgs, *istlComm_);
418 amg_ =
new AMG(fineOperator, coarsenCriterion, smootherArgs);
422 const Simulator &simulator_;
427 OwnerOverlapCopyCommunication *istlComm_;
429 OverlappingMatrix *overlappingMatrix_;
430 OverlappingVector *overlappingb_;
431 OverlappingVector *overlappingx_;
Uses communication on the grid to find the initial seed list of indices.
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.
Provides a linear solver backend using the parallel algebraic multi-grid (AMG) linear solver from DUN...
Definition: parallelamgbackend.hh:60
Base class for all convergence criteria which only defines an virtual API.
Definition: convergencecriterion.hh:51
static void registerParameters()
Definition: parallelamgbackend.hh:145
An overlap aware ISTL scalar product.
SET_INT_PROP(NumericModel, GridGlobalRefinements, 0)
GridManager & gridManager()
Return a reference to the grid manager of simulation.
Definition: simulator.hh:152
Convergence criterion which looks at the weighted absolute value of the residual. ...
Definition: weightedresidreductioncriterion.hh:54
~ParallelAmgBackend()
Definition: parallelamgbackend.hh:142
Definition: cartesianindexmapper.hh:31
NEW_PROP_TAG(Grid)
The type of the DUNE grid.
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
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:64
void setStructureMatrix(const Matrix &M)
Set the structure of the linear system of equations to be solved.
Definition: parallelamgbackend.hh:162
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:176
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:623
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
ParallelAmgBackend(const Simulator &simulator)
Definition: parallelamgbackend.hh:132
Provides the magic behind the eWoms property system.
An overlap aware block vector.
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:229
Implements a generic linear solver abstraction.
An overlap aware preconditioner for any ISTL linear solver.
bool solve(const Matrix &M, Vector &x, Vector &b)
Actually solve the linear system of equations.
Definition: parallelamgbackend.hh:173
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:95