5 #include <dune/common/unused.hh> 6 #include <dune/common/version.hh> 8 #include <dune/istl/novlpschwarz.hh> 9 #include <dune/istl/schwarz.hh> 10 #include <dune/istl/paamg/amg.hh> 11 #include <dune/istl/paamg/fastamg.hh> 12 #include <dune/istl/paamg/kamg.hh> 14 #include <amdis/Initfile.hpp> 15 #include <amdis/common/ConceptsBase.hpp> 16 #include <amdis/linearalgebra/istl/Communication.hpp> 17 #include <amdis/linearalgebra/istl/ISTLPreconCreator.hpp> 18 #include <amdis/linearalgebra/istl/precompiled/Preconditioners.hpp> 19 #include <amdis/linearalgebra/istl/precompiled/Solvers.hpp> 26 template <
template <
class...>
class AMGSolver,
class Traits>
29 template <
class Smoother>
30 using SmootherArgs =
typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
34 template <
class Traits>
37 using PrecBase =
typename Traits::Prec;
39 template <
class Smoother,
class LinOp,
class Criterion,
class Comm>
40 static std::unique_ptr<PrecBase>
41 create([[maybe_unused]] std::string prefix, LinOp
const& linOp, Criterion
const& criterion, SmootherArgs<Smoother>
const& smootherArgs, Comm
const& comm)
43 using Solver = Dune::Amg::AMG<LinOp, typename Traits::X, Smoother, Comm>;
44 return std::make_unique<Solver>(linOp, criterion, smootherArgs, comm);
54 template <
class Traits>
57 using PrecBase =
typename Traits::Prec;
59 template <
class Smoother,
class LinOp,
class Criterion,
class SmootherArgs,
class Comm>
60 static std::unique_ptr<PrecBase>
61 create(std::string
const& prefix, LinOp
const& linOp, Criterion
const& criterion,
62 SmootherArgs
const& smootherArgs, Comm
const& comm)
64 return createImpl1<Smoother>(prefix, linOp, criterion, smootherArgs, comm, Dune::PriorityTag<5>{});
69 template <
class Smoother,
class LinOp,
class Criterion>
70 static std::unique_ptr<PrecBase>
71 createImpl1(std::string prefix, LinOp
const& linOp, Criterion
const& criterion,
72 [[maybe_unused]] SmootherArgs<Smoother>
const& smootherArgs,
73 Dune::Amg::SequentialInformation
const& comm, Dune::PriorityTag<2>)
75 bool symmetric = Parameters::get<bool>(prefix +
"->symmetric").value_or(
true);
77 using Solver = Dune::Amg::FastAMG<LinOp, typename Traits::X, Dune::Amg::SequentialInformation>;
78 return std::make_unique<Solver>(linOp, criterion, criterion, symmetric, comm);
81 template <
class Smoother,
class LinOp,
class Criterion,
class SmootherArgs,
class Comm>
82 static std::unique_ptr<PrecBase>
83 createImpl1(std::string, LinOp
const&, Criterion
const&, SmootherArgs
const&, Comm
const&,
86 error_exit(
"Currently only sequential FastAMG supported.");
102 template <
class Traits>
105 using X =
typename Traits::X;
106 using PrecBase =
typename Traits::Prec;
108 template <
class Smoother,
class LinOp,
class Criterion,
class Comm>
109 static std::unique_ptr<PrecBase>
110 create(std::string prefix, LinOp
const& linOp, Criterion
const& criterion,
111 [[maybe_unused]] SmootherArgs<Smoother>
const& smootherArgs, Comm
const& comm)
113 std::string solver = Parameters::get<std::string>(prefix +
"->krylov solver").value_or(
"default");
115 std::size_t maxLevelKrylovSteps = 3;
116 Parameters::get(prefix +
"->krylov solver->maxLevelKrylovSteps", maxLevelKrylovSteps);
117 double minDefectReduction = 1.e-1;
118 Parameters::get(prefix +
"->krylov solver->minDefectReduction", minDefectReduction);
120 if (solver ==
"pcg" || solver ==
"default")
122 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::GeneralizedPCGSolver<X>>;
123 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
125 else if (solver ==
"bicgstab" || solver ==
"bcgs")
127 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::BiCGSTABSolver<X>>;
128 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
130 else if (solver ==
"minres")
132 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::MINRESSolver<X>>;
133 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
138 else if (solver ==
"gmres")
140 using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::RestartedGMResSolver<X>>;
141 return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
146 error_exit(
"Unknown coarse space solver {} given for parameter `{}`.", solver, prefix +
"->coarse space solver");
175 template <
template <
class...>
class AMGSolver,
class Traits>
179 using M =
typename Traits::M;
180 using X =
typename Traits::X;
181 using Y =
typename Traits::Y;
183 using Interface = Dune::Preconditioner<X,Y>;
184 using LinOperator =
typename Traits::LinOpCreator::Interface;
186 using SolverCategory =
typename Dune::SolverCategory::Category;
188 static const bool is_arithmetic = std::is_arithmetic_v<typename Dune::FieldTraits<M>::field_type>;
189 using Norm = std::conditional_t<is_arithmetic, Dune::Amg::FirstDiagonal, Dune::Amg::RowSum>;
191 using SymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<M,Norm>>;
192 using UnSymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<M,Norm>>;
196 void init(std::string
const& prefix)
override 199 initParameters(prefix);
203 std::unique_ptr<Interface>
204 create(M
const& mat,
typename Traits::Comm
const& comm)
const override 206 return createImpl0(Dune::SolverCategory::category(comm), mat, comm.get());
211 template <
class Comm>
212 std::unique_ptr<Interface>
213 createImpl0(SolverCategory cat, M
const& mat, Comm
const& comm)
const 215 if (smoother_ ==
"diag" || smoother_ ==
"jacobi" || smoother_ ==
"default") {
216 return createImpl1<Dune::SeqJac<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
221 else if (smoother_ ==
"sor") {
222 return createImpl1<Dune::SeqSOR<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
227 else if (smoother_ ==
"gs" || smoother_ ==
"gauss_seidel") {
228 return createImpl1<Dune::SeqGS<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
233 else if (smoother_ ==
"richardson") {
234 return createImpl1<Dune::Richardson<X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
237 else if (smoother_ ==
"ssor") {
238 return createImpl1<Dune::SeqSSOR<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
241 error_exit(
"Unknown smoother type '{}' given for parameter '{}'", smoother_, prefix_ +
"->smoother");
248 template <
class Smoother>
249 std::unique_ptr<Interface>
250 createImpl1([[maybe_unused]] SolverCategory cat, M
const& mat,
251 Dune::Amg::SequentialInformation
const& comm, Dune::PriorityTag<2>)
const 253 assert(cat == SolverCategory::sequential);
254 using LinOp = Dune::MatrixAdapter<M,X,Y>;
255 LinOp* linOpPtr =
new LinOp(mat);
256 linOperator_.reset(linOpPtr);
257 return createImpl2<Smoother>(*linOpPtr, comm);
260 template <
class Smoother,
class Comm>
261 std::unique_ptr<Interface>
262 createImpl1(SolverCategory cat, M
const& mat, Comm
const& comm, Dune::PriorityTag<1>)
const 265 case SolverCategory::sequential:
267 return createImpl1<Smoother>(cat, mat, Dune::Amg::SequentialInformation{}, Dune::PriorityTag<5>{});
270 case SolverCategory::overlapping:
272 using LinOp = Dune::OverlappingSchwarzOperator<M,X,Y,Comm>;
273 using ParSmoother = Dune::BlockPreconditioner<X,Y,Comm,Smoother>;
274 LinOp* linOpPtr =
new LinOp(mat, comm);
275 linOperator_.reset(linOpPtr);
276 return createImpl2<ParSmoother>(*linOpPtr, comm);
278 case SolverCategory::nonoverlapping:
280 using LinOp = Dune::NonoverlappingSchwarzOperator<M,X,Y,Comm>;
281 using ParSmoother = Dune::NonoverlappingBlockPreconditioner<Comm,Smoother>;
282 LinOp* linOpPtr =
new LinOp(mat, comm);
283 linOperator_.reset(linOpPtr);
284 return createImpl2<ParSmoother>(*linOpPtr, comm);
288 error_exit(
"Invalid solver category for AMGSolver\n");
295 template <
class Smoother,
class LinOp,
class Comm>
296 std::unique_ptr<Interface> createImpl2(LinOp
const& linOp, Comm
const& comm)
const 298 SmootherArgs<Smoother> smootherArgs;
299 initSmootherParameters(prefix_ +
"->smoother", smootherArgs);
301 return symmetricAggregation_
302 ? createImpl3<Smoother>(linOp, SymCriterion(parameters_), smootherArgs, comm)
303 : createImpl3<Smoother>(linOp, UnSymCriterion(parameters_), smootherArgs, comm);
308 template <
class Smoother,
class LinOp,
class Criterion,
class Comm>
309 std::unique_ptr<Interface>
310 createImpl3(LinOp
const& linOp, Criterion
const& criterion,
311 SmootherArgs<Smoother>
const& smootherArgs, Comm
const& comm)
const 314 return Creator::template create<Smoother>(prefix_,linOp,criterion,smootherArgs,comm);
318 void initParameters(std::string
const& prefix)
321 smoother_ =
"default";
325 parameters_.setDebugLevel(Parameters::get<int>(prefix +
"->debugLevel").value_or(2));
327 parameters_.setNoPreSmoothSteps(Parameters::get<std::size_t>(prefix +
"->preSmoothSteps").value_or(2));
329 parameters_.setNoPostSmoothSteps(Parameters::get<std::size_t>(prefix +
"->postSmoothSteps").value_or(2));
331 parameters_.setGamma(Parameters::get<std::size_t>(prefix +
"->gamma").value_or(1));
333 parameters_.setAdditive(Parameters::get<bool>(prefix +
"->additive").value_or(
false));
335 symmetricAggregation_ = Parameters::get<bool>(prefix +
"->symmetric aggregation").value_or(
true);
337 initCoarseningParameters(prefix +
"->coarsening");
338 initAggregationParameters(prefix +
"->aggregation");
339 initDependencyParameters(prefix +
"->dependency");
342 void initCoarseningParameters(std::string
const& prefix)
345 parameters_.setMaxLevel(Parameters::get<int>(prefix +
"->maxLevel").value_or(100));
347 parameters_.setCoarsenTarget(Parameters::get<int>(prefix +
"->coarsenTarget").value_or(1000));
349 parameters_.setMinCoarsenRate(Parameters::get<double>(prefix +
"->minCoarsenRate").value_or(1.2));
351 parameters_.setProlongationDampingFactor(Parameters::get<double>(prefix +
"->dampingFactor").value_or(1.6));
354 void initAggregationParameters(std::string
const& prefix)
357 parameters_.setMaxDistance(Parameters::get<std::size_t>(prefix +
"->maxDistance").value_or(2));
359 parameters_.setMinAggregateSize(Parameters::get<std::size_t>(prefix +
"->minAggregateSize").value_or(4));
361 parameters_.setMaxAggregateSize(Parameters::get<std::size_t>(prefix +
"->maxAggregateSize").value_or(6));
363 parameters_.setMaxConnectivity(Parameters::get<std::size_t>(prefix +
"->maxConnectivity").value_or(15));
365 parameters_.setSkipIsolated(Parameters::get<bool>(prefix +
"->skipIsolated").value_or(
false));
368 void initDependencyParameters(std::string
const& prefix)
371 parameters_.setAlpha(Parameters::get<double>(prefix +
"->alpha").value_or(1.0/3.0));
373 parameters_.setBeta(Parameters::get<double>(prefix +
"->beta").value_or(1.e-5));
377 void initSmootherParameters(std::string
const& prefix, SA& smootherArgs)
const 380 Parameters::get(prefix +
"->relaxationFactor", smootherArgs.relaxationFactor);
385 std::string smoother_;
386 Dune::Amg::Parameters parameters_;
387 bool symmetricAggregation_ =
true;
389 mutable std::shared_ptr<LinOperator> linOperator_;
Definition: AdaptiveGrid.hpp:373
Definition: AdaptBase.hpp:6
std::unique_ptr< Interface > create(M const &mat, typename Traits::Comm const &comm) const override
Implements ISTLPreconCreatorBase::create.
Definition: AMGPrecon.hpp:204
Definition: AMGPrecon.hpp:27
Base class for precon creators,.
Definition: ISTLPreconCreator.hpp:39
static std::optional< T > get(std::string const &key)
Get parameter-values from parameter-tree.
Definition: Initfile.hpp:25
void init(std::string const &prefix) override
Prepare the preconditioner for the creation.
Definition: AMGPrecon.hpp:196
Definition: AMGPrecon.hpp:176