27 #ifndef EWOMS_PARALLEL_AMG_BACKEND_HH
28 #define EWOMS_PARALLEL_AMG_BACKEND_HH
34 #include <dune/istl/paamg/amg.hh>
35 #include <dune/istl/paamg/pinfo.hh>
36 #include <dune/istl/owneroverlapcopy.hh>
42 template <
class TypeTag>
46 namespace Properties {
54 SET_INT_PROP(ParallelAmgLinearSolver, AmgCoarsenTarget, 5000);
69 template <
class TypeTag>
70 class ParallelAmgBackend :
public ParallelBaseBackend<TypeTag>
72 typedef ParallelBaseBackend<TypeTag> ParentType;
74 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
75 typedef typename GET_PROP_TYPE(TypeTag, LinearSolverScalar) LinearSolverScalar;
77 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
78 typedef typename GET_PROP_TYPE(TypeTag, Overlap) Overlap;
80 typedef typename ParentType::ParallelOperator ParallelOperator;
81 typedef typename ParentType::OverlappingVector OverlappingVector;
82 typedef typename ParentType::ParallelPreconditioner ParallelPreconditioner;
83 typedef typename ParentType::ParallelScalarProduct ParallelScalarProduct;
86 typedef Dune::FieldVector<LinearSolverScalar, numEq> VectorBlock;
87 typedef Dune::FieldMatrix<LinearSolverScalar, numEq, numEq> MatrixBlock;
89 typedef Dune::BCRSMatrix<MatrixBlock> Matrix;
90 typedef Dune::BlockVector<VectorBlock> Vector;
94 typedef Dune::SeqSOR<Matrix, Vector, Vector> SequentialSmoother;
101 typedef Dune::OwnerOverlapCopyCommunication<Ewoms::Linear::Index>
102 OwnerOverlapCopyCommunication;
103 typedef Dune::OverlappingSchwarzOperator<Matrix,
106 OwnerOverlapCopyCommunication> FineOperator;
107 typedef Dune::OverlappingSchwarzScalarProduct<Vector,
108 OwnerOverlapCopyCommunication> FineScalarProduct;
109 typedef Dune::BlockPreconditioner<Vector,
111 OwnerOverlapCopyCommunication,
112 SequentialSmoother> ParallelSmoother;
113 typedef Dune::Amg::AMG<FineOperator,
116 OwnerOverlapCopyCommunication> AMG;
118 typedef Dune::MatrixAdapter<Matrix, Vector, Vector> FineOperator;
119 typedef Dune::SeqScalarProduct<Vector> FineScalarProduct;
120 typedef SequentialSmoother ParallelSmoother;
121 typedef Dune::Amg::AMG<FineOperator, Vector, ParallelSmoother> AMG;
124 typedef BiCGStabSolver<ParallelOperator,
126 AMG> RawLinearSolver;
129 ParallelAmgBackend(
const Simulator& simulator)
130 : ParentType(simulator)
133 static void registerParameters()
138 "The maximum residual error which the linear solver tolerates"
139 " without giving up");
141 "The coarsening target for the agglomerations of "
142 "the AMG preconditioner");
148 std::shared_ptr<AMG> preparePreconditioner_()
153 istlComm_ = std::make_shared<OwnerOverlapCopyCommunication>(MPI_COMM_WORLD);
154 setupAmgIndexSet_(this->overlappingMatrix_->overlap(), istlComm_->indexSet());
155 istlComm_->remoteIndices().template rebuild<false>();
160 fineOperator_ = std::make_shared<FineOperator>(*this->overlappingMatrix_, *istlComm_);
162 fineOperator_ = std::make_shared<FineOperator>(*this->overlappingMatrix_);
170 void cleanupPreconditioner_()
173 std::shared_ptr<RawLinearSolver> prepareSolver_(ParallelOperator& parOperator,
174 ParallelScalarProduct& parScalarProduct,
177 const auto& gridView = this->simulator_.
gridView();
178 typedef CombinedCriterion<OverlappingVector, decltype(gridView.comm())> CCC;
180 Scalar linearSolverTolerance =
EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance);
181 Scalar linearSolverAbsTolerance = this->simulator_.
model().newtonMethod().tolerance() / 10.0;
183 convCrit_.reset(
new CCC(gridView.comm(),
184 linearSolverTolerance,
185 linearSolverAbsTolerance,
188 auto bicgstabSolver =
189 std::make_shared<RawLinearSolver>(parPreCond, *convCrit_, parScalarProduct);
192 if (parOperator.overlap().myRank() == 0)
194 bicgstabSolver->setVerbosity(verbosity);
195 bicgstabSolver->setMaxIterations(
EWOMS_GET_PARAM(TypeTag,
int, LinearSolverMaxIterations));
196 bicgstabSolver->setLinearOperator(&parOperator);
197 bicgstabSolver->setRhs(this->overlappingb_);
199 return bicgstabSolver;
202 bool runSolver_(std::shared_ptr<RawLinearSolver> solver)
203 {
return solver->apply(*this->overlappingx_); }
205 void cleanupSolver_()
209 template <
class ParallelIndexSet>
210 void setupAmgIndexSet_(
const Overlap& overlap, ParallelIndexSet& istlIndices)
212 typedef Dune::OwnerOverlapCopyAttributeSet GridAttributes;
213 typedef Dune::OwnerOverlapCopyAttributeSet::AttributeSet GridAttributeSet;
216 istlIndices.beginResize();
217 for (Index curIdx = 0;
static_cast<size_t>(curIdx) < overlap.numDomestic(); ++curIdx) {
218 GridAttributeSet gridFlag =
219 overlap.iAmMasterOf(curIdx)
220 ? GridAttributes::owner
221 : GridAttributes::copy;
225 bool isShared = overlap.isInOverlap(curIdx);
227 assert(curIdx == overlap.globalToDomestic(overlap.domesticToGlobal(curIdx)));
228 istlIndices.add(overlap.domesticToGlobal(curIdx),
229 Dune::ParallelLocalIndex<GridAttributeSet>(
static_cast<size_t>(curIdx),
233 istlIndices.endResize();
243 if (this->simulator_.
gridManager().gridView().comm().rank() == 0)
246 typedef typename Dune::Amg::SmootherTraits<ParallelSmoother>::Arguments SmootherArgs;
248 SmootherArgs smootherArgs;
249 smootherArgs.iterations = 1;
250 smootherArgs.relaxationFactor = 1.0;
258 CoarsenCriterion<Dune::Amg::SymmetricCriterion<Matrix, Dune::Amg::FrobeniusNorm> >
261 CoarsenCriterion coarsenCriterion(15, coarsenTarget);
262 coarsenCriterion.setDefaultValuesAnisotropic(GridView::dimension,
265 coarsenCriterion.setDebugLevel(1);
267 coarsenCriterion.setDebugLevel(0);
270 coarsenCriterion.setMinCoarsenRate(1.05);
272 coarsenCriterion.setAccumulate(Dune::Amg::atOnceAccu);
273 coarsenCriterion.setSkipIsolated(
false);
277 amg_ = std::make_shared<AMG>(*fineOperator_, coarsenCriterion, smootherArgs, *istlComm_);
279 amg_ = std::make_shared<AMG>(*fineOperator_, coarsenCriterion, smootherArgs);
283 std::unique_ptr<ConvergenceCriterion<OverlappingVector> > convCrit_;
285 std::shared_ptr<FineOperator> fineOperator_;
286 std::shared_ptr<AMG> amg_;
289 std::shared_ptr<OwnerOverlapCopyCommunication> istlComm_;
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
Provides a linear solver backend using the parallel algebraic multi-grid (AMG) linear solver from DUN...
Definition: parallelamgbackend.hh:43
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
#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
#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
Implements a preconditioned stabilized BiCG 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
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
const GridView & gridView() const
Return the grid view for which the simulation is done.
Definition: simulator.hh:183
#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.