1#ifndef DUNE_COPASI_SOLVER_ISTL_FACTORY_ITERATIVE_HH
2#define DUNE_COPASI_SOLVER_ISTL_FACTORY_ITERATIVE_HH
9#include <dune/istl/blocklevel.hh>
10#include <dune/istl/preconditioner.hh>
11#include <dune/istl/solver.hh>
12#include <dune/istl/solvers.hh>
14#include <dune/common/classname.hh>
15#include <dune/common/exceptions.hh>
17#include <fmt/ranges.h>
19#define DUNE_COPASI_DEFAULT_ITERATIVE_SOLVER "BiCGSTAB"
24template<
class X,
class Y,
class A>
26 std::shared_ptr<IterativeSolver<X, Y>>(
const std::shared_ptr<LinearOperator<X, Y>>&,
27 const std::shared_ptr<ScalarProduct<X>>&,
28 const std::shared_ptr<Preconditioner<X, Y>>&,
31template<
class X,
class Y,
class A>
36template<
class Solver,
class Alloc>
38defaultIterativeSolverFactory()
40 using X =
typename Solver::domain_type;
41 using Y =
typename Solver::range_type;
42 return [](
const std::shared_ptr<LinearOperator<X, Y>>& op,
43 const std::shared_ptr<ScalarProduct<X>>& sp,
44 const std::shared_ptr<Preconditioner<X, Y>>& prec,
45 const ParameterTree& config,
46 const Alloc& alloc) -> std::shared_ptr<IterativeSolver<X, Y>> {
48 using SolverAlloc =
typename std::allocator_traits<Alloc>::template rebind_alloc<Solver>;
49 SolverAlloc salloc(alloc);
52 const auto& conv_cond = config.sub(
"convergence_condition");
53 auto reduction = conv_cond.get(
"relative_tolerance", 1e-8);
54 auto [min_it, max_it] = conv_cond.get(
"iteration_range", std::array<std::size_t, 2>{ 1, 500 });
57 "Minimum number of iterations '{}' cannot be set for solver '{}. It has been set to '1'.'",
60 auto verbosity = config.get(
"verbosity", 0);
63 if constexpr (std::derived_from<Solver, RestartedGMResSolver<X, Y>>) {
64 auto restart = config.get(
"restart", 40);
65 return std::allocate_shared<Solver>(salloc, op, sp, prec, reduction, restart, max_it, verbosity);
66 }
else if constexpr (std::derived_from<Solver, GeneralizedPCGSolver<X>> or std::derived_from<Solver, RestartedFCGSolver<X>>) {
67 auto restart = config.get(
"restart", 10);
68 return std::allocate_shared<Solver>(salloc, op, sp, prec, reduction, max_it, verbosity, restart);
70 return std::allocate_shared<Solver>(salloc, op, sp, prec, reduction, max_it, verbosity);
76template<Concept::LinearOperator O,
class Alloc>
77const IterativeSolverRegistry<typename O::domain_type, typename O::range_type, Alloc>&
80 using X =
typename O::domain_type;
81 using Y =
typename O::range_type;
85 static std::once_flag flag;
86 std::call_once(flag, [] {
88 Impl::defaultIterativeSolverFactory<LoopSolver<X>, Alloc>());
89 registry.
define(
"Gradient",
90 Impl::defaultIterativeSolverFactory<GradientSolver<X>, Alloc>());
92 Impl::defaultIterativeSolverFactory<CGSolver<X>, Alloc>());
93 registry.
define(
"BiCGSTAB",
94 Impl::defaultIterativeSolverFactory<BiCGSTABSolver<X>, Alloc>());
96 Impl::defaultIterativeSolverFactory<MINRESSolver<X>, Alloc>());
97 registry.
define(
"RestartedGMRes",
98 Impl::defaultIterativeSolverFactory<RestartedGMResSolver<X, Y>, Alloc>());
99 registry.
define(
"RestartedFlexibleGMRes",
100 Impl::defaultIterativeSolverFactory<RestartedFlexibleGMResSolver<X, Y>, Alloc>());
101 registry.
define(
"GeneralizedPCG",
102 Impl::defaultIterativeSolverFactory<GeneralizedPCGSolver<X>, Alloc>());
103 registry.
define(
"RestartedFCG",
104 Impl::defaultIterativeSolverFactory<RestartedFCGSolver<X>, Alloc>());
105 registry.
define(
"CompleteFCG",
106 Impl::defaultIterativeSolverFactory<CompleteFCGSolver<X>, Alloc>());
111template<Concept::LinearOperator Op,
class Alloc = std::allocator<Op>>
112std::shared_ptr<IterativeSolver<typename Op::domain_type, typename Op::range_type>>
114 const std::shared_ptr<Op>& op,
115 const std::shared_ptr<ScalarProduct<typename Op::domain_type>>& sp,
116 const std::shared_ptr<Preconditioner<typename Op::domain_type, typename Op::range_type>>& prec,
117 const ParameterTree& config,
118 const Alloc& alloc = Alloc())
123 auto& registry = getIterativeSolverRegistry<Op, Alloc>();
125 if (config.get(
"verbosity", 1) > 0)
126 spdlog::info(
"Creating iterative solver with type '{}'", type_name);
127 if (registry.contains(type_name))
128 return registry.create(type_name, op, sp, prec, config, alloc);
132 "The key '{}' is not a known iterative solver type. Allowed types are {}",
simple wrapper that captures and exposes the keys of Dune::ParameterizedObjectFactory
Definition: parameterized_object.hh:14
void define(Key const &key, T &&t)
Definition: parameterized_object.hh:21
#define DUNE_COPASI_DEFAULT_ITERATIVE_SOLVER
Definition: iterative.hh:19
Definition: block_jacobi.hh:14
const IterativeSolverRegistry< typename O::domain_type, typename O::range_type, Alloc > & getIterativeSolverRegistry()
Definition: iterative.hh:78
std::shared_ptr< IterativeSolver< typename Op::domain_type, typename Op::range_type > > makeIterativeSolver(const std::shared_ptr< Op > &op, const std::shared_ptr< ScalarProduct< typename Op::domain_type > > &sp, const std::shared_ptr< Preconditioner< typename Op::domain_type, typename Op::range_type > > &prec, const ParameterTree &config, const Alloc &alloc=Alloc())
Definition: iterative.hh:113
std::shared_ptr< IterativeSolver< X, Y > >(const std::shared_ptr< LinearOperator< X, Y > > &, const std::shared_ptr< ScalarProduct< X > > &, const std::shared_ptr< Preconditioner< X, Y > > &, const ParameterTree &, const A &) IterativeSolverFactorySignature
Definition: iterative.hh:30
auto format_exception(Exception &&e, fmt::format_string< Args... > format, Args &&... args)
Definition: exceptions.hh:23