97 using OverlappingVector =
typename ParentType::OverlappingVector;
102 using VectorBlock = Dune::FieldVector<LinearSolverScalar, numEq>;
103 using MatrixBlock =
typename SparseMatrixAdapter::MatrixBlock;
104 using IstlMatrix =
typename SparseMatrixAdapter::IstlMatrix;
106 using Vector = Dune::BlockVector<VectorBlock>;
110 using SequentialSmoother = Dune::SeqSOR<IstlMatrix, Vector, Vector>;
116 using OwnerOverlapCopyCommunication = Dune::OwnerOverlapCopyCommunication<Opm::Linear::Index>;
117 using FineOperator = Dune::OverlappingSchwarzOperator<IstlMatrix,
120 OwnerOverlapCopyCommunication>;
121 using FineScalarProduct = Dune::OverlappingSchwarzScalarProduct<Vector,
122 OwnerOverlapCopyCommunication>;
123 using ParallelSmoother = Dune::BlockPreconditioner<Vector,
125 OwnerOverlapCopyCommunication,
127 using AMG = Dune::Amg::AMG<FineOperator,
130 OwnerOverlapCopyCommunication>;
132 using FineOperator = Dune::MatrixAdapter<IstlMatrix, Vector, Vector>;
133 using FineScalarProduct = Dune::SeqScalarProduct<Vector>;
134 using ParallelSmoother = SequentialSmoother;
135 using AMG = Dune::Amg::AMG<FineOperator, Vector, ParallelSmoother>;
142 static_assert(std::is_same<SparseMatrixAdapter, IstlSparseMatrixAdapter<MatrixBlock> >::value,
143 "The ParallelAmgBackend linear solver backend requires the IstlSparseMatrixAdapter");
150 static void registerParameters()
154 Parameters::Register<Parameters::LinearSolverMaxError<Scalar>>
155 (
"The maximum residual error which the linear solver tolerates "
156 "without giving up");
157 Parameters::Register<Parameters::AmgCoarsenTarget>
158 (
"The coarsening target for the agglomerations of "
159 "the AMG preconditioner");
165 std::shared_ptr<AMG> preparePreconditioner_()
170 istlComm_ = std::make_shared<OwnerOverlapCopyCommunication>(MPI_COMM_WORLD);
171 setupAmgIndexSet_(this->overlappingMatrix_->overlap(), istlComm_->indexSet());
172 istlComm_->remoteIndices().template rebuild<false>();
177 fineOperator_ = std::make_shared<FineOperator>(*this->overlappingMatrix_, *istlComm_);
179 fineOperator_ = std::make_shared<FineOperator>(*this->overlappingMatrix_);
187 void cleanupPreconditioner_()
190 std::shared_ptr<RawLinearSolver> prepareSolver_(ParallelOperator& parOperator,
191 ParallelScalarProduct& parScalarProduct,
194 const auto& gridView = this->simulator_.gridView();
197 Scalar linearSolverTolerance = Parameters::Get<Parameters::LinearSolverTolerance<Scalar>>();
198 Scalar linearSolverAbsTolerance = Parameters::Get<Parameters::LinearSolverAbsTolerance<Scalar>>();
199 if (linearSolverAbsTolerance < 0.0) {
200 linearSolverAbsTolerance = this->simulator_.model().newtonMethod().tolerance()/100.0;
203 convCrit_.reset(
new CCC(gridView.comm(),
204 linearSolverTolerance,
205 linearSolverAbsTolerance,
208 auto bicgstabSolver =
209 std::make_shared<RawLinearSolver>(parPreCond, *convCrit_, parScalarProduct);
212 if (parOperator.overlap().myRank() == 0)
213 verbosity = Parameters::Get<Parameters::LinearSolverVerbosity>();
214 bicgstabSolver->setVerbosity(verbosity);
215 bicgstabSolver->setMaxIterations(Parameters::Get<Parameters::LinearSolverMaxIterations>());
216 bicgstabSolver->setLinearOperator(&parOperator);
217 bicgstabSolver->setRhs(this->overlappingb_);
219 return bicgstabSolver;
222 std::pair<bool,int> runSolver_(std::shared_ptr<RawLinearSolver> solver)
224 bool converged = solver->apply(*this->overlappingx_);
225 return std::make_pair(converged,
int(solver->report().iterations()));
228 void cleanupSolver_()
232 template <
class ParallelIndexSet>
233 void setupAmgIndexSet_(
const Overlap& overlap, ParallelIndexSet& istlIndices)
235 using GridAttributes = Dune::OwnerOverlapCopyAttributeSet;
236 using GridAttributeSet = Dune::OwnerOverlapCopyAttributeSet::AttributeSet;
239 istlIndices.beginResize();
240 for (
Index curIdx = 0;
static_cast<size_t>(curIdx) < overlap.numDomestic(); ++curIdx) {
241 GridAttributeSet gridFlag =
242 overlap.iAmMasterOf(curIdx)
243 ? GridAttributes::owner
244 : GridAttributes::copy;
248 bool isShared = overlap.isInOverlap(curIdx);
250 assert(curIdx == overlap.globalToDomestic(overlap.domesticToGlobal(curIdx)));
251 istlIndices.add(overlap.domesticToGlobal(curIdx),
252 Dune::ParallelLocalIndex<GridAttributeSet>(
static_cast<size_t>(curIdx),
256 istlIndices.endResize();
261 template <
typename C>
262 auto setUseFixedOrder(C criterion,
bool booleanValue) ->
decltype(criterion.setUseFixedOrder(booleanValue))
264 return criterion.setUseFixedOrder(booleanValue);
266 template <
typename C>
267 void setUseFixedOrder(C, ...)
278 if (this->simulator_.vanguard().gridView().comm().rank() == 0)
279 verbosity = Parameters::Get<Parameters::LinearSolverVerbosity>();
281 using SmootherArgs =
typename Dune::Amg::SmootherTraits<ParallelSmoother>::Arguments;
283 SmootherArgs smootherArgs;
284 smootherArgs.iterations = 1;
285 smootherArgs.relaxationFactor = 1.0;
292 using CoarsenCriterion = Dune::Amg::
293 CoarsenCriterion<Dune::Amg::SymmetricCriterion<IstlMatrix, Dune::Amg::FrobeniusNorm> >;
294 int coarsenTarget = Parameters::Get<Parameters::AmgCoarsenTarget>();
295 CoarsenCriterion coarsenCriterion(15, coarsenTarget);
296 coarsenCriterion.setDefaultValuesAnisotropic(GridView::dimension,
299 coarsenCriterion.setDebugLevel(1);
301 coarsenCriterion.setDebugLevel(0);
304 coarsenCriterion.setMinCoarsenRate(1.05);
306 coarsenCriterion.setAccumulate(Dune::Amg::atOnceAccu);
307 coarsenCriterion.setSkipIsolated(
false);
308 setUseFixedOrder(coarsenCriterion,
true);
312 amg_ = std::make_shared<AMG>(*fineOperator_, coarsenCriterion, smootherArgs, *istlComm_);
314 amg_ = std::make_shared<AMG>(*fineOperator_, coarsenCriterion, smootherArgs);
318 std::unique_ptr<ConvergenceCriterion<OverlappingVector> > convCrit_;
320 std::shared_ptr<FineOperator> fineOperator_;
321 std::shared_ptr<AMG> amg_;
324 std::shared_ptr<OwnerOverlapCopyCommunication> istlComm_;