1#ifndef DUNE_FEM_SOLVER_ISTL_HH
2#define DUNE_FEM_SOLVER_ISTL_HH
19#include <dune/common/exceptions.hh>
20#include <dune/common/ftraits.hh>
21#include <dune/common/hybridutilities.hh>
22#include <dune/common/typelist.hh>
23#include <dune/common/typeutilities.hh>
26#include <dune/geometry/dimension.hh>
29#include <dune/istl/bcrsmatrix.hh>
30#include <dune/istl/operators.hh>
31#include <dune/istl/paamg/amg.hh>
32#include <dune/istl/preconditioners.hh>
33#include <dune/istl/schwarz.hh>
34#include <dune/istl/solvers.hh>
59 template<
class M,
class X,
class Y>
60 struct ConstructionTraits<SeqILDL<M, X, Y>>
62 using Arguments = DefaultConstructionArgs<SeqILDL<M, X, Y>>;
64 static inline auto construct (Arguments& args) -> std::shared_ptr<SeqILDL<M, X, Y>>
66 return std::make_shared<SeqILDL<M, X, Y>>(args.getMatrix(), args.getArgs().relaxationFactor);
80 template<
class DiscreteFunctionSpace>
81 class HierarchicalDiscreteFunction;
83 template<
class DomainFunction,
class RangeFunction>
84 struct ISTLLinearOperator;
94 template<
class DiscreteFunction >
95 using VectorType = std::decay_t<decltype(std::declval<const DiscreteFunction&>().dofVector().array())>;
102 template<
class LinearOperator>
105 using Type = std::decay_t<decltype(std::declval<const LinearOperator&>().exportMatrix())>;
108 template<
class DomainFunction,
class RangeFunction>
109 struct __MatrixType<
Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>>
111 using Type = Dune::BCRSMatrix<typename Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>::LittleBlockType>;
115 template<
class LinearOperator>
116 using MatrixType =
typename __MatrixType<LinearOperator>::Type;
122 template<
template<
class,
class,
class,
int...>
class Prec,
class Op,
int ... l>
123 using __FillPrecondType = Prec<
typename Op::matrix_type,
typename Op::domain_type,
typename Op::range_type, l ...>;
130 struct IsBCRSMatrix : std::false_type {};
132 template<
class B,
class A>
133 struct IsBCRSMatrix<
Dune::BCRSMatrix<B, A>> : std::true_type {};
140 enum Symmetry :
bool { unsymmetric =
false, symmetric =
true };
147 template<
class T,
class =
void>
148 struct UniformLeafLevelType;
151 using UniformLeafLevel =
typename UniformLeafLevelType<T>::Type;
153 template<
class V1,
class... V>
154 struct UniformLeafLevelType<
Dune::MultiTypeBlockVector<V1, V...>, std::enable_if_t<Std::are_all_same<UniformLeafLevel<V1>, UniformLeafLevel<V>...>::value>>
156 using Type = std::integral_constant<int, UniformLeafLevel<V1>::value + 1>;
159 template<
class B,
class A>
160 struct UniformLeafLevelType<
Dune::BlockVector<B, A>>
162 using Type = std::integral_constant<int, UniformLeafLevel<B>::value + 1>;
165 template<
class K,
int n>
166 struct UniformLeafLevelType<
Dune::FieldVector<K, n>>
168 using Type = std::integral_constant<int, 0>;
171 template<
class R1,
class... R>
172 struct UniformLeafLevelType<
Dune::MultiTypeBlockMatrix< R1, R...>, std::enable_if_t<Std::are_all_same<UniformLeafLevel<R1>, UniformLeafLevel<R>...>::value>>
175 using Type = std::integral_constant<int, UniformLeafLevel<R1>::value>;
178 template<
class B,
class A>
179 struct UniformLeafLevelType<
Dune::BCRSMatrix<B, A>>
181 using Type = std::integral_constant<int, UniformLeafLevel<B>::value + 1>;
184 template<
class K,
int m,
int n>
185 struct UniformLeafLevelType<
Dune::FieldMatrix<K, m, n>>
187 using Type = std::integral_constant<int, 0>;
195 template<
class M,
class X,
class Y,
196 std::enable_if_t<IsBCRSMatrix<M>::value && (UniformLeafLevel<M>::value == 1),
int> = 0>
197 inline static decltype(
auto) namedSmootherTypes (PriorityTag<2>)
199 return std::make_tuple(std::make_pair(std::string(
"jacobi"), Dune::MetaType<Dune::SeqJac< M, X, Y, 1>>{}),
200 std::make_pair(std::string(
"gauss-seidel"), Dune::MetaType<Dune::SeqGS< M, X, Y, 1>>{}),
201 std::make_pair(std::string(
"sor"), Dune::MetaType<Dune::SeqSOR<M, X, Y, 1>>{}),
202 std::make_pair(std::string(
"ssor"), Dune::MetaType<Dune::SeqSSOR<M, X, Y, 1>>{}),
203 std::make_pair(std::string(
"ilu"), Dune::MetaType<Dune::SeqILU< M, X, Y>>{}),
204 std::make_pair(std::string(
"ildl"), Dune::MetaType<Dune::SeqILDL<M, X, Y>>{}));
207 template<
class M,
class X,
class Y,
208 std::enable_if_t<(UniformLeafLevel<M>::value >= 0),
int> = 0>
209 inline static decltype(
auto) namedSmootherTypes (PriorityTag<1>)
211 return std::make_tuple(std::make_pair(std::string(
"jacobi"), Dune::MetaType<Dune::SeqJac<M, X, Y, UniformLeafLevel<M>::value>>{}),
212 std::make_pair(std::string(
"gauss-seidel"), Dune::MetaType<Dune::SeqGS<M, X, Y, UniformLeafLevel<M>::value>>{}),
213 std::make_pair(std::string(
"sor"), Dune::MetaType<Dune::SeqSOR<M, X, Y, UniformLeafLevel<M>::value>>{}),
214 std::make_pair(std::string(
"ssor"), Dune::MetaType<Dune::SeqSSOR<M, X, Y, UniformLeafLevel<M>::value>>{}));
217 template<
class M,
class X,
class Y>
218 inline static decltype(
auto) namedSmootherTypes (PriorityTag<0>)
220 return std::make_tuple();
223 template<
class M,
class X,
class Y>
224 inline static decltype(
auto) namedSmootherTypes ()
226 return namedSmootherTypes<M, X, Y>(PriorityTag<42>{});
230 inline static decltype(
auto) namedAMGNormTypes ()
232 return std::make_tuple(std::make_pair(std::string(
"firstdiagonal"), Dune::MetaType<Dune::Amg::FirstDiagonal>{}),
233 std::make_pair(std::string(
"rowsum"), Dune::MetaType<Dune::Amg::RowSum>{}),
234 std::make_pair(std::string(
"frobenius"), Dune::MetaType<Dune::Amg::FrobeniusNorm>{}),
235 std::make_pair(std::string(
"one"), Dune::MetaType<Dune::Amg::AlwaysOneNorm>{}));
243 class SolverParameter :
public LocalParameter<Fem::SolverParameter, SolverParameter>
245 using BaseType = LocalParameter<Fem::SolverParameter, SolverParameter>;
248 using BaseType::parameter;
249 using BaseType::keyPrefix;
252 template<
class... T,
class F>
253 void getEnum (
const std::string& key,
254 const std::tuple<std::pair<std::string, T>...>& choices,
255 const std::string& defaultValue,
258 const std::string& value = parameter().getValue(key, defaultValue);
259 bool success =
false;
261 if (value != choice.first)
273 Hybrid::forEach(choices, [&values] (
const auto& choice) { values += (values.empty() ?
"'" :
", '") + choice.first +
"'"; });
274 DUNE_THROW(ParameterInvalid,
"Parameter '" << key <<
"' invalid (choices are: " << values <<
").");
278 using BaseType::gmresRestart;
279 using BaseType::relaxation;
280 using BaseType::preconditionerIteration;
281 using BaseType::preconditionerLevel;
282 using BaseType::preconditionMethod;
283 using BaseType::preconditionMethodTable;
284 using BaseType::verbose;
285 using BaseType::setVerbose;
289 : BaseType(
"istl.", parameter)
293 : BaseType(keyPrefix, parameter)
296 SolverParameter (
const Fem::SolverParameter& parameter)
297 : SolverParameter(parameter.keyPrefix(), parameter.parameter())
300 virtual int errorMeasure ()
const override
302 if (errorMeasure_ < 0)
303 errorMeasure_ = BaseType::errorMeasure();
304 return errorMeasure_;
307 virtual int verbosity ()
const
309 const std::string verbosityTable[] = {
"off",
"on",
"full" };
310 return parameter().getEnum(keyPrefix() +
"verbosity", verbosityTable, 0);
313 virtual int restart ()
const {
return parameter().getValue(keyPrefix() +
"restart", 20); }
314 virtual int fcgMmax ()
const {
return parameter().getValue(keyPrefix() +
"fcg.mmax", 10); }
316 virtual double precRelaxation ()
const {
return parameter().getValue(keyPrefix() +
"preconditioner.relax", 1.0); }
317 virtual int precIterations ()
const {
return parameter().getValue(keyPrefix() +
"preconditioner.iterations", 1); }
319 virtual int precType ()
const
321 const std::string types[] = {
"richardson",
"smoother",
"amg" };
322 return parameter().getEnum(keyPrefix() +
"preconditioner.type", types, 1);
325 template<
class... T,
class F>
326 void precSmoother (
const std::tuple<std::pair<std::string, T>...>& choices,
const std::string& defaultChoice, F&& f)
const
328 getEnum(keyPrefix() +
"preconditioner.smoother", choices, defaultChoice, std::forward<F>(f));
331 virtual int iluFillin ()
const {
return parameter().getValue(keyPrefix() +
"ilu.fillin", 0); }
332 virtual bool iluReorder ()
const {
return parameter().getValue(keyPrefix() +
"ilu.reorder",
false); }
334 virtual int amgMaxLevel ()
const {
return parameter().getValue(keyPrefix() +
"amg.maxlevel", 100); }
335 virtual int amgCoarsenTarget ()
const {
return parameter().getValue(keyPrefix() +
"amg.coarsentarget", 1000); }
336 virtual double amgMinCoarsenRate ()
const {
return parameter().getValue(keyPrefix() +
"amg.mincoarsenrate", 1.2); }
337 virtual double amgProlongDamping ()
const {
return parameter().getValue(keyPrefix() +
"amg.prolongation.dampingfactor", 1.6); }
338 virtual int amgDebugLevel ()
const {
return parameter().getValue(keyPrefix() +
"amg.debuglevel", 0); }
339 virtual int amgPreSmoothSteps ()
const {
return parameter().getValue(keyPrefix() +
"amg.presmoothsteps", 2); }
340 virtual int amgPostSmoothSteps ()
const {
return parameter().getValue(keyPrefix() +
"amg.postsmoothsteps", 2); }
341 virtual bool amgAdditive ()
const {
return parameter().getValue(keyPrefix() +
"amg.additive",
false); }
342 virtual double amgAlpha ()
const {
return parameter().getValue(keyPrefix() +
"amg.alpha", 1.0/3.0); }
343 virtual double amgBeta ()
const {
return parameter().getValue(keyPrefix() +
"amg.beta", 1.0e-5); }
345 virtual int amgCycle ()
const
347 const std::string cycles[] = {
"v-cycle",
"w-cycle" };
348 return 1 + parameter().getEnum(keyPrefix() +
"amg.cycle", cycles, 0);
351 virtual int amgAccumulate ()
const
353 const std::string modes[] = {
"none",
"once",
"successive" };
354 return parameter().getEnum(keyPrefix() +
"amg.accumulate", modes, 2);
357 virtual std::size_t amgAggregationDimension ()
const {
return parameter().getValue<std::size_t>(keyPrefix() +
"amg.aggregation.dimension", 2); }
358 virtual std::size_t amgAggregationDistance ()
const {
return parameter().getValue<std::size_t>(keyPrefix() +
"amg.aggregation.distance", 2); }
359 virtual std::size_t amgAggregationMinSize ()
const {
return parameter().getValue<std::size_t>(keyPrefix() +
"amg.aggregation.min", 4); }
360 virtual std::size_t amgAggregationMaxSize ()
const {
return parameter().getValue<std::size_t>(keyPrefix() +
"amg.aggregation.max", 6); }
361 virtual std::size_t amgAggregationConnectivity ()
const {
return parameter().getValue<std::size_t>(keyPrefix() +
"amg.aggregation.connectivity", 15); }
362 virtual bool amgAggregationSkipIsolated ()
const {
return parameter().getValue(keyPrefix() +
"amg.aggregation.skipisolated",
false); }
364 virtual int amgAggregationMode ()
const
366 const std::string modes[] = {
"isotropic",
"anisotropic",
"parameter" };
367 return parameter().getEnum(keyPrefix() +
"amg.aggregation.mode", modes, 0);
370 template<
class... T,
class F>
371 void amgNorm (
const std::tuple<std::pair<std::string, T>...>& choices,
const std::string& defaultChoice, F&& f)
const
373 getEnum(keyPrefix() +
"amg.norm", choices, defaultChoice, std::forward<F>(f));
377 mutable int errorMeasure_ = -1;
385 template<
class AssembledOperator,
class Communication,
386 std::enable_if_t<SupportsAMG<Communication>::value,
int> = 0>
387 inline auto makeAMGPreconditioner (
const std::shared_ptr<AssembledOperator> op,
const Communication& comm, Symmetry symmetry,
388 const SolverParameter& parameter = {})
389 -> std::shared_ptr<Dune::Preconditioner<typename AssembledOperator::domain_type, typename AssembledOperator::range_type>>
391 using matrix_type =
typename AssembledOperator::matrix_type;
392 using domain_type =
typename AssembledOperator::domain_type;
393 using range_type =
typename AssembledOperator::range_type;
394 using real_type = real_t<typename AssembledOperator::field_type>;
396 std::shared_ptr<Dune::Preconditioner<domain_type, range_type>> preconditioner;
397 const auto smoothers = namedSmootherTypes<matrix_type, domain_type, range_type>();
398 parameter.precSmoother(smoothers,
"jacobi", [op, &comm, symmetry, ¶meter, &preconditioner] (
auto type) {
399 using SmootherType = Dune::BlockPreconditioner<domain_type, range_type, Communication,
typename decltype(type)::type>;
400 using AMG = Dune::Amg::AMG<AssembledOperator, domain_type, SmootherType, Communication>;
402 Dune::Amg::DefaultSmootherArgs <real_type> smootherArgs;
403 smootherArgs.relaxationFactor = parameter.precRelaxation();
404 smootherArgs.iterations = parameter.precIterations();
406 Dune::Amg::Parameters amgParams(
407 parameter.amgMaxLevel(),
408 parameter.amgCoarsenTarget(),
409 parameter.amgMinCoarsenRate(),
410 parameter.amgProlongDamping(),
411 static_cast<Dune::Amg::AccumulationMode
>(parameter.amgAccumulate())
415 switch(parameter.amgAggregationMode())
418 amgParams.setDefaultValuesIsotropic(parameter.amgAggregationDimension(), parameter.amgAggregationDistance());
422 amgParams.setDefaultValuesAnisotropic(parameter.amgAggregationDimension(), parameter.amgAggregationDistance());
426 amgParams.setMaxDistance(parameter.amgAggregationDistance());
427 amgParams.setMinAggregateSize(parameter.amgAggregationMinSize());
428 amgParams.setMaxAggregateSize(parameter.amgAggregationMaxSize());
432 amgParams.setMaxConnectivity(parameter.amgAggregationConnectivity());
433 amgParams.setSkipIsolated(parameter.amgAggregationSkipIsolated());
435 amgParams.setDebugLevel(parameter.amgDebugLevel());
436 amgParams.setNoPreSmoothSteps(parameter.amgPreSmoothSteps());
437 amgParams.setNoPostSmoothSteps(parameter.amgPostSmoothSteps());
438 amgParams.setAlpha(parameter.amgAlpha());
439 amgParams.setBeta(parameter.amgBeta());
440 amgParams.setGamma(parameter.amgCycle());
441 amgParams.setAdditive(parameter.amgAdditive());
443 parameter.amgNorm(namedAMGNormTypes(),
"rowsum", [op, &comm, symmetry, &amgParams, &smootherArgs, &preconditioner](
auto type){
444 if (symmetry == symmetric)
446 Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion< matrix_type,
typename decltype(type)::type>> criterion(amgParams);
447 preconditioner.reset(
new AMG(*op, criterion, smootherArgs, comm), [op] (AMG* p) {
delete p; });
451 Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<matrix_type,
typename decltype(type)::type>> criterion( amgParams );
452 preconditioner.reset(
new AMG(*op, criterion, smootherArgs, comm), [op] (AMG* p) {
delete p; });
456 return preconditioner;
459 template<
class AssembledOperator,
class Communication,
460 std::enable_if_t<!SupportsAMG<Communication>::value,
int> = 0>
461 inline auto makeAMGPreconditioner (
const std::shared_ptr<AssembledOperator> op,
const Communication& comm, Symmetry symmetry,
462 const SolverParameter& parameter = {})
463 -> std::shared_ptr<Dune::Preconditioner<typename AssembledOperator::domain_type, typename AssembledOperator::range_type>>
465 DUNE_THROW(Dune::InvalidStateException,
"Communication does not support AMG.");
474 inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
475 std::shared_ptr<Dune::Richardson<typename Op::domain_type, typename Op::range_type>>& preconditioner,
476 const SolverParameter& parameter = {})
478 using Preconditioner = Dune::Richardson<typename Op::domain_type, typename Op::range_type>;
479 preconditioner.reset(
new Preconditioner(parameter.precRelaxation()),
480 [op] ( Preconditioner* p) { delete p; });
483 template<
class Op,
int l>
484 inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
485 std::shared_ptr<__FillPrecondType<Dune::SeqJac, Op, l>>& preconditioner,
486 const SolverParameter& parameter = {})
488 using Preconditioner = __FillPrecondType<Dune::SeqJac, Op, l>;
489 preconditioner.reset(
new Preconditioner(op->getmat(), parameter.precIterations(), parameter.precRelaxation()),
490 [op] (Preconditioner* p) { delete p; });
493 template<
class Op,
int l>
494 inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
495 std::shared_ptr<__FillPrecondType<Dune::SeqSOR, Op, l>>& preconditioner,
496 const SolverParameter& parameter = {})
498 using Preconditioner = __FillPrecondType<Dune::SeqSOR, Op, l>;
499 preconditioner.reset(
new Preconditioner(op->getmat(), parameter.precIterations(), parameter.precRelaxation()),
500 [op] (Preconditioner* p) { delete p; });
503 template<
class Op,
int l>
504 inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
505 std::shared_ptr<__FillPrecondType<Dune::SeqSSOR, Op, l>>& preconditioner,
506 const SolverParameter& parameter = {})
508 using Preconditioner = __FillPrecondType<Dune::SeqSSOR, Op, l>;
509 preconditioner.reset(
new Preconditioner(op->getmat(), parameter.precIterations(), parameter.precRelaxation()),
510 [op] (Preconditioner* p) { delete p; });
514 inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
515 std::shared_ptr<__FillPrecondType<Dune::SeqILU, Op>>& preconditioner,
516 const SolverParameter& parameter = {})
518 using Preconditioner = __FillPrecondType<Dune::SeqILU, Op>;
519 preconditioner.reset(
new Preconditioner(op->getmat(), parameter.iluFillin(), parameter.precRelaxation(), parameter.iluReorder()),
520 [op] (Preconditioner* p) { delete p; });
524 inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
525 std::shared_ptr<__FillPrecondType<Dune::SeqILDL, Op>>& preconditioner,
526 const SolverParameter& parameter = {})
528 using Preconditioner = __FillPrecondType<Dune::SeqILDL, Op>;
529 preconditioner.reset(
new Preconditioner(op->getmat(), parameter.precRelaxation()),
530 [op] (Preconditioner* p) { delete p; });
538 template<
class SeqPreconditioner,
class Communication>
539 inline auto makeBlockPreconditioner(std::shared_ptr<SeqPreconditioner> seqPreconditioner,
const Communication& comm)
540 -> std::shared_ptr<Dune::Preconditioner<typename SeqPreconditioner::domain_type, typename SeqPreconditioner::range_type>>
542 using domain_type =
typename SeqPreconditioner::domain_type;
543 using range_type =
typename SeqPreconditioner::range_type;
545 using BlockPreconditioner = Dune::BlockPreconditioner<domain_type, range_type, Communication, SeqPreconditioner>;
547 return std::make_shared<BlockPreconditioner>(seqPreconditioner, comm);
555 template<
class Op,
class Communication>
556 inline auto makePreconditioner (std::shared_ptr<Op> op,
const Communication& comm, Symmetry symmetry,
557 const SolverParameter& parameter = {})
558 -> std::shared_ptr<Dune::Preconditioner<typename Op::domain_type, typename Op::range_type>>
560 using domain_type =
typename Op::domain_type;
561 using range_type =
typename Op::range_type;
563 const auto smoothers = namedSmootherTypes<typename Op::matrix_type, domain_type, range_type>();
564 std::shared_ptr<Dune::Preconditioner<domain_type, range_type>> preconditioner;
566 switch (parameter.precType())
570 std::shared_ptr<Dune::Richardson<domain_type, range_type>> seqPreconditioner;
571 makeSequentialPreconditioner(op, seqPreconditioner, parameter);
572 preconditioner = makeBlockPreconditioner(seqPreconditioner, comm);
577 parameter.precSmoother(smoothers,
"jacobi", [op, &comm, ¶meter, &preconditioner] (
auto type) {
578 std::shared_ptr<
typename decltype(type)::type> seqPreconditioner;
579 makeSequentialPreconditioner(op, seqPreconditioner, parameter);
580 preconditioner = makeBlockPreconditioner(seqPreconditioner, comm);
585 preconditioner = makeAMGPreconditioner(op, comm, symmetry, parameter);
589 DUNE_THROW(Dune::InvalidStateException,
"Invalid ISTL preconditioner type selected.");
591 return preconditioner;
600 class InverseOperator final
601 :
public Dune::Fem::Operator<typename LinearOperator::RangeFunctionType, typename LinearOperator::DomainFunctionType>
603 static_assert(std::is_same< typename LinearOperator::DomainFunctionType, typename LinearOperator::RangeFunctionType>::value,
604 "Domain function and range function must have the same type.");
607 using LinearOperatorType = LinearOperator;
608 using OperatorType = LinearOperatorType;
610 using DiscreteFunctionType =
typename LinearOperatorType::DomainFunctionType;
611 using DiscreteFunctionSpaceType =
typename DiscreteFunctionType::DiscreteFunctionSpaceType;
613 using RealType = real_t<typename DiscreteFunctionType::RangeFieldType>;
615 using SolverParameterType = SolverParameter;
618 using vector_type = VectorType<DiscreteFunctionType>;
619 using matrix_type = MatrixType<LinearOperatorType>;
622 using CommunicationType = Communication<DiscreteFunctionSpaceType>;
624 using AssembledLinearOperatorType = Dune::OverlappingSchwarzOperator<matrix_type, vector_type, vector_type, CommunicationType>;
625 using PreconditionerType = Dune::Preconditioner<vector_type, vector_type>;
626 using ScalarProductType = Dune::ScalarProduct<vector_type>;
628 using PreconditionerFactory = std::function<std::shared_ptr<PreconditionerType>(std::shared_ptr<AssembledLinearOperatorType>,
const CommunicationType&, Symmetry,
const SolverParameterType&)>;
630 InverseOperator (
const SolverParameterType& parameter = {})
631 : InverseOperator(makePreconditioner<AssembledLinearOperatorType, CommunicationType>, parameter)
634 InverseOperator (PreconditionerFactory factory,
const SolverParameterType& parameter = {})
635 : preconditionerFactory_{std::move(factory)},
636 parameter_{std::make_shared<SolverParameterType>(parameter)}
639 InverseOperator (
const LinearOperatorType& op,
const SolverParameterType& parameter = {})
640 : InverseOperator(parameter)
645 InverseOperator (
const LinearOperatorType& op, PreconditionerFactory factory,
const SolverParameterType& parameter = {})
646 : InverseOperator(std::move(factory), parameter)
651 void operator() (
const DiscreteFunctionType& u, DiscreteFunctionType& w)
const override
654 vector_type b(u.dofVector().array());
655 if (parameter_->errorMeasure() == 0)
657 vector_type residuum(b);
658 linearOperator_->applyscaleadd(-1.0, w.dofVector().array(), residuum);
659 RealType res = scalarProduct_->norm(residuum);
660 RealType reduction = (res > 0) ? parameter_->tolerance() / res : 1e-3;
661 solver_->apply(w.dofVector().array(), b, reduction, result_);
664 solver_->apply(w.dofVector().array(), b, result_);
667 void bind (
const LinearOperatorType& op)
669 buildCommunication(op.domainSpace(), Dune::SolverCategory::overlapping, communication_);
671 linearOperator_.reset(
new AssembledLinearOperatorType(op.exportMatrix(), *communication_));
672 scalarProduct_.reset(
new Dune::OverlappingSchwarzScalarProduct<vector_type, CommunicationType>(*communication_));
675 preconditioner_ = preconditionerFactory_(linearOperator_, *communication_, symmetry, parameter());
678 auto reduction = parameter().tolerance();
679 auto maxIterations = parameter().maxIterations();
680 auto verbosity = op.domainSpace().gridPart().comm().rank() == 0 ? parameter().verbosity() : 0;
682 switch (parameter().solverMethod({SolverParameterType::cg,
683 SolverParameterType::bicgstab,
684 SolverParameterType::gmres,
685 SolverParameterType::minres,
686 SolverParameterType::gradient,
687 SolverParameterType::loop},
688 {
"pcg",
"fcg",
"fgmres"},
689 (symmetry == symmetric ? 0 : 1)))
692 solver_.reset(
new Dune::RestartedFlexibleGMResSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, parameter().restart(), maxIterations, verbosity));
696 solver_.reset(
new Dune::RestartedFCGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity, parameter().fcgMmax()));
700 solver_.reset(
new Dune::GeneralizedPCGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity, parameter().restart()));
704 solver_.reset(
new Dune::CGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
708 solver_.reset(
new Dune::BiCGSTABSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
712 solver_.reset(
new Dune::RestartedGMResSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, parameter().restart(), maxIterations, verbosity));
716 solver_.reset(
new Dune::MINRESSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
720 solver_.reset(
new Dune::GradientSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
724 solver_.reset(
new Dune::LoopSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
732 preconditioner_.reset();
733 scalarProduct_.reset();
734 linearOperator_.reset();
735 communication_.reset();
738 auto parameter () -> SolverParameterType& {
return *parameter_; }
739 void setParamters (
const SolverParameterType& parameter) { parameter_ = std::make_shared<SolverParameterType>(parameter); }
741 int iterations ()
const {
return result_.iterations; }
742 bool converged ()
const {
return result_.converged; }
744 void setMaxIterations (
int maxIterations) { parameter().setMaxIterations(maxIterations); }
747 PreconditionerFactory preconditionerFactory_;
748 std::shared_ptr<SolverParameterType> parameter_;
750 std::shared_ptr<CommunicationType> communication_;
751 std::shared_ptr<AssembledLinearOperatorType> linearOperator_;
752 std::shared_ptr<ScalarProductType> scalarProduct_;
753 std::shared_ptr<PreconditionerType> preconditioner_;
754 std::shared_ptr<Dune::InverseOperator<vector_type, vector_type>> solver_;
756 mutable Dune::InverseOperatorResult result_;
765 class OperatorPreconditionerFactory
767 static_assert(std::is_same<typename Operator::DomainFunctionType, typename Operator::RangeFunctionType>::value,
768 "Domain function and range function must have the same type.");
770 using LinearOperatorType =
typename Operator::JacobianOperatorType;
771 using DiscreteFunctionType =
typename LinearOperatorType::DomainFunctionType;
773 using vector_type = VectorType<DiscreteFunctionType>;
774 using matrix_type = MatrixType<LinearOperatorType>;
777 using DiscreteFunctionSpaceType =
typename DiscreteFunctionType::DiscreteFunctionSpaceType;
779 using CommunicationType = Communication<DiscreteFunctionSpaceType>;
781 using PreconditionerType = Dune::Preconditioner<vector_type, vector_type>;
783 template<
class... Args>
784 OperatorPreconditionerFactory (
const DiscreteFunctionSpaceType& space, Args&& ...args)
785 : op_(
std::forward<Args>(args)...),
786 zero_(
"zero", space),
787 jacobian_(
"preconditioner", space, space)
790 template<
class AssembledOperator>
791 auto operator() (std::shared_ptr<AssembledOperator>,
const CommunicationType& comm, Symmetry,
792 const SolverParameter& parameter = {})
const
793 -> std::shared_ptr<PreconditionerType>
795 using AssembledLinearOperatorType = Dune::OverlappingSchwarzOperator<matrix_type, vector_type, vector_type, CommunicationType>;
796 op_.jacobian(zero_, jacobian_);
797 return makePreconditioner(std::make_shared<AssembledLinearOperatorType>(jacobian_.matrix(), comm), comm, symmetry, parameter);
800 auto op () const -> const Operator& {
return op_; }
801 auto op () -> Operator& {
return op_; }
805 DiscreteFunctionType zero_;
806 mutable LinearOperatorType jacobian_;
Definition: bindguard.hh:11
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
Dune::OwnerOverlapCopyCommunication< std::size_t, typename DiscreteFunctionSpace::BlockMapperType::GlobalKeyType > OwnerOverlapCopyCommunication
Definition: owneroverlapcopy.hh:34
void buildCommunication(const DiscreteFunctionSpace &dfSpace, Dune::SolverCategory::Category solverCategory, std::shared_ptr< FemCommunication< DiscreteFunctionSpace > > &communication)
Definition: fem.hh:143
static ParameterContainer & container()
Definition: io/parameter.hh:193
abstract operator
Definition: operator.hh:34