DUNE-FEM (unstable)

istl.hh
1#ifndef DUNE_FEM_SOLVER_ISTL_HH
2#define DUNE_FEM_SOLVER_ISTL_HH
3
4#if HAVE_DUNE_ISTL
5
6#include <cassert>
7
8// stl includes
9#include <limits>
10#include <map>
11#include <memory>
12#include <string>
13#include <type_traits>
14#include <tuple>
15#include <utility>
16#include <vector>
17
18// dune-common includes
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>
24
25// dune-geometry includes
26#include <dune/geometry/dimension.hh>
27
28// dune-istl includes
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>
35
36// dune-fem includes
37#include <dune/fem/common/utility.hh>
38#include <dune/fem/io/parameter.hh>
39#include <dune/fem/operator/common/operator.hh>
40#include <dune/fem/solver/parameter.hh>
41#include <dune/fem/space/mapper/parallel.hh>
42
43// local includes
44#include <dune/fem/solver/communication/fem.hh>
45#include <dune/fem/solver/communication/hierarchical.hh>
46#include <dune/fem/solver/communication/owneroverlapcopy.hh>
47
48
49
50namespace Dune
51{
52
53 namespace Amg
54 {
55
56 // ConstructionTraits for SeqILDL
57 // ------------------------------
58
59 template<class M, class X, class Y>
60 struct ConstructionTraits<SeqILDL<M, X, Y>>
61 {
62 using Arguments = DefaultConstructionArgs<SeqILDL<M, X, Y>>;
63
64 static inline auto construct (Arguments& args) -> std::shared_ptr<SeqILDL<M, X, Y>>
65 {
66 return std::make_shared<SeqILDL<M, X, Y>>(args.getMatrix(), args.getArgs().relaxationFactor);
67 }
68 };
69
70 } // namespace Amg
71
72
73
74 namespace Fem
75 {
76
77 // External Forward Declarations
78 // -----------------------------
79
80 template<class DiscreteFunctionSpace>
81 class HierarchicalDiscreteFunction;
82
83 template<class DomainFunction, class RangeFunction>
84 struct ISTLLinearOperator;
85
86
87
88 namespace ISTL
89 {
90
91 // VectorType
92 // ----------
93
94 template< class DiscreteFunction >
95 using VectorType = std::decay_t<decltype(std::declval<const DiscreteFunction&>().dofVector().array())>;
96
97
98
99 // MatrixType
100 // ----------
101
102 template<class LinearOperator>
103 struct __MatrixType
104 {
105 using Type = std::decay_t<decltype(std::declval<const LinearOperator&>().exportMatrix())>;
106 };
107
108 template<class DomainFunction, class RangeFunction>
109 struct __MatrixType<Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>>
110 {
111 using Type = Dune::BCRSMatrix<typename Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>::LittleBlockType>;
112 };
113
114
115 template<class LinearOperator>
116 using MatrixType = typename __MatrixType<LinearOperator>::Type;
117
118
119 // __FillPrecondType
120 // -----------------
121
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 ...>;
124
125
126 // IsBCRSMatrix
127 // ------------
128
129 template<class M>
130 struct IsBCRSMatrix : std::false_type {};
131
132 template<class B, class A>
133 struct IsBCRSMatrix<Dune::BCRSMatrix<B, A>> : std::true_type {};
134
135
136
137 // Symmetry
138 // --------
139
140 enum Symmetry : bool { unsymmetric = false, symmetric = true };
141
142
143
144 // UniformLeafLevel
145 // ----------------
146
147 template<class T, class = void>
148 struct UniformLeafLevelType;
149
150 template<class T>
151 using UniformLeafLevel = typename UniformLeafLevelType<T>::Type;
152
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>>
155 {
156 using Type = std::integral_constant<int, UniformLeafLevel<V1>::value + 1>;
157 };
158
159 template<class B, class A>
160 struct UniformLeafLevelType<Dune::BlockVector<B, A>>
161 {
162 using Type = std::integral_constant<int, UniformLeafLevel<B>::value + 1>;
163 };
164
165 template<class K, int n>
166 struct UniformLeafLevelType<Dune::FieldVector<K, n>>
167 {
168 using Type = std::integral_constant<int, 0>;
169 };
170
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>>
173 {
174 // Note: The rows of MultiTypeBlockMatrix are of type MultiTypeBlockVector, which already increases the level
175 using Type = std::integral_constant<int, UniformLeafLevel<R1>::value>;
176 };
177
178 template<class B, class A>
179 struct UniformLeafLevelType<Dune::BCRSMatrix<B, A>>
180 {
181 using Type = std::integral_constant<int, UniformLeafLevel<B>::value + 1>;
182 };
183
184 template<class K, int m, int n>
185 struct UniformLeafLevelType<Dune::FieldMatrix<K, m, n>>
186 {
187 using Type = std::integral_constant<int, 0>;
188 };
189
190
191
192 // namedSmootherTypes
193 // ------------------
194
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>)
198 {
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>>{}));
205 }
206
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>)
210 {
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>>{}));
215 }
216
217 template<class M, class X, class Y>
218 inline static decltype(auto) namedSmootherTypes (PriorityTag<0>)
219 {
220 return std::make_tuple();
221 }
222
223 template<class M, class X, class Y>
224 inline static decltype(auto) namedSmootherTypes ()
225 {
226 return namedSmootherTypes<M, X, Y>(PriorityTag<42>{});
227 }
228
229
230 inline static decltype(auto) namedAMGNormTypes ()
231 {
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>{}));
236 }
237
238
239
240 // SolverParameter
241 // ---------------
242
243 class SolverParameter : public LocalParameter<Fem::SolverParameter, SolverParameter>
244 {
245 using BaseType = LocalParameter<Fem::SolverParameter, SolverParameter>;
246
247 public:
248 using BaseType::parameter;
249 using BaseType::keyPrefix;
250
251 private:
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,
256 F&& f) const
257 {
258 const std::string& value = parameter().getValue(key, defaultValue);
259 bool success = false;
260 Hybrid::forEach(choices, [&f, &value, &success] (const auto& choice) {
261 if (value != choice.first)
262 return;
263
264 assert(!success);
265 success = true;
266
267 f(choice.second);
268 });
269 if (success)
270 return;
271
272 std::string values;
273 Hybrid::forEach(choices, [&values] (const auto& choice) { values += (values.empty() ? "'" : ", '") + choice.first + "'"; });
274 DUNE_THROW(ParameterInvalid, "Parameter '" << key << "' invalid (choices are: " << values << ").");
275 }
276
277 // hide functions
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;
286
287 public:
288 SolverParameter (const ParameterReader& parameter = Parameter::container())
289 : BaseType("istl.", parameter)
290 {}
291
292 SolverParameter (const std::string& keyPrefix, const ParameterReader& parameter = Parameter::container())
293 : BaseType(keyPrefix, parameter)
294 {}
295
296 SolverParameter (const Fem::SolverParameter& parameter)
297 : SolverParameter(parameter.keyPrefix(), parameter.parameter())
298 {}
299
300 virtual int errorMeasure () const override
301 {
302 if (errorMeasure_ < 0)
303 errorMeasure_ = BaseType::errorMeasure();
304 return errorMeasure_;
305 }
306
307 virtual int verbosity () const
308 {
309 const std::string verbosityTable[] = { "off", "on", "full" };
310 return parameter().getEnum(keyPrefix() + "verbosity", verbosityTable, 0);
311 }
312
313 virtual int restart () const { return parameter().getValue(keyPrefix() + "restart", 20); }
314 virtual int fcgMmax () const { return parameter().getValue(keyPrefix() + "fcg.mmax", 10); }
315
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); }
318
319 virtual int precType () const
320 {
321 const std::string types[] = { "richardson", "smoother", "amg" };
322 return parameter().getEnum(keyPrefix() + "preconditioner.type", types, 1);
323 }
324
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
327 {
328 getEnum(keyPrefix() + "preconditioner.smoother", choices, defaultChoice, std::forward<F>(f));
329 }
330
331 virtual int iluFillin () const { return parameter().getValue(keyPrefix() + "ilu.fillin", 0); }
332 virtual bool iluReorder () const { return parameter().getValue(keyPrefix() + "ilu.reorder", false); }
333
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); }
344
345 virtual int amgCycle () const
346 {
347 const std::string cycles[] = { "v-cycle", "w-cycle" };
348 return 1 + parameter().getEnum(keyPrefix() + "amg.cycle", cycles, 0);
349 }
350
351 virtual int amgAccumulate () const
352 {
353 const std::string modes[] = { "none", "once", "successive" };
354 return parameter().getEnum(keyPrefix() + "amg.accumulate", modes, 2);
355 }
356
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); }
363
364 virtual int amgAggregationMode () const
365 {
366 const std::string modes[] = { "isotropic", "anisotropic", "parameter" };
367 return parameter().getEnum(keyPrefix() + "amg.aggregation.mode", modes, 0);
368 }
369
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
372 {
373 getEnum(keyPrefix() + "amg.norm", choices, defaultChoice, std::forward<F>(f));
374 }
375
376 private:
377 mutable int errorMeasure_ = -1;
378 };
379
380
381
382 // makeAMGPreconditioner
383 // ---------------------
384
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>>
390 {
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>;
395
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, &parameter, &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>;
401
402 Dune::Amg::DefaultSmootherArgs <real_type> smootherArgs;
403 smootherArgs.relaxationFactor = parameter.precRelaxation();
404 smootherArgs.iterations = parameter.precIterations();
405
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())
412 );
413
414 // parameters
415 switch(parameter.amgAggregationMode())
416 {
417 case 0:
418 amgParams.setDefaultValuesIsotropic(parameter.amgAggregationDimension(), parameter.amgAggregationDistance()); // ...
419 break;
420
421 case 1:
422 amgParams.setDefaultValuesAnisotropic(parameter.amgAggregationDimension(), parameter.amgAggregationDistance()); // ...
423 break;
424
425 case 2:
426 amgParams.setMaxDistance(parameter.amgAggregationDistance());
427 amgParams.setMinAggregateSize(parameter.amgAggregationMinSize());
428 amgParams.setMaxAggregateSize(parameter.amgAggregationMaxSize());
429 break;
430 }
431
432 amgParams.setMaxConnectivity(parameter.amgAggregationConnectivity());
433 amgParams.setSkipIsolated(parameter.amgAggregationSkipIsolated());
434
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());
442
443 parameter.amgNorm(namedAMGNormTypes(), "rowsum", [op, &comm, symmetry, &amgParams, &smootherArgs, &preconditioner](auto type){
444 if (symmetry == symmetric)
445 {
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; });
448 }
449 else
450 {
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; });
453 }
454 });
455 });
456 return preconditioner;
457 }
458
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>>
464 {
465 DUNE_THROW(Dune::InvalidStateException, "Communication does not support AMG.");
466 }
467
468
469
470 // makeSequentialPreconditioner
471 // ----------------------------
472
473 template<class Op>
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 = {})
477 {
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; });
481 }
482
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 = {})
487 {
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; });
491 }
492
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 = {})
497 {
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; });
501 }
502
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 = {})
507 {
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; });
511 }
512
513 template<class Op>
514 inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
515 std::shared_ptr<__FillPrecondType<Dune::SeqILU, Op>>& preconditioner,
516 const SolverParameter& parameter = {})
517 {
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; });
521 }
522
523 template<class Op>
524 inline void makeSequentialPreconditioner (std::shared_ptr<Op> op,
525 std::shared_ptr<__FillPrecondType<Dune::SeqILDL, Op>>& preconditioner,
526 const SolverParameter& parameter = {})
527 {
528 using Preconditioner = __FillPrecondType<Dune::SeqILDL, Op>;
529 preconditioner.reset(new Preconditioner(op->getmat(), parameter.precRelaxation()),
530 [op] (Preconditioner* p) { delete p; });
531 }
532
533
534
535 // makeBlockPreconditioner
536 // -----------------------
537
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>>
541 {
542 using domain_type = typename SeqPreconditioner::domain_type;
543 using range_type = typename SeqPreconditioner::range_type;
544
545 using BlockPreconditioner = Dune::BlockPreconditioner<domain_type, range_type, Communication, SeqPreconditioner>;
546
547 return std::make_shared<BlockPreconditioner>(seqPreconditioner, comm);
548 }
549
550
551
552 // makePreconditioner
553 // ------------------
554
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>>
559 {
560 using domain_type = typename Op::domain_type;
561 using range_type = typename Op::range_type;
562
563 const auto smoothers = namedSmootherTypes<typename Op::matrix_type, domain_type, range_type>();
564 std::shared_ptr<Dune::Preconditioner<domain_type, range_type>> preconditioner;
565
566 switch (parameter.precType())
567 {
568 case 0:
569 {
570 std::shared_ptr<Dune::Richardson<domain_type, range_type>> seqPreconditioner;
571 makeSequentialPreconditioner(op, seqPreconditioner, parameter);
572 preconditioner = makeBlockPreconditioner(seqPreconditioner, comm);
573 }
574 break;
575
576 case 1:
577 parameter.precSmoother(smoothers, "jacobi", [op, &comm, &parameter, &preconditioner] (auto type) {
578 std::shared_ptr<typename decltype(type)::type> seqPreconditioner;
579 makeSequentialPreconditioner(op, seqPreconditioner, parameter);
580 preconditioner = makeBlockPreconditioner(seqPreconditioner, comm);
581 } );
582 break;
583
584 case 2:
585 preconditioner = makeAMGPreconditioner(op, comm, symmetry, parameter);
586 break;
587
588 default:
589 DUNE_THROW(Dune::InvalidStateException, "Invalid ISTL preconditioner type selected.");
590 }
591 return preconditioner;
592 }
593
594
595
596 // InverseOperator
597 // ---------------
598
599 template<class LinearOperator, Symmetry symmetry = symmetric, template<class> class Communication = OwnerOverlapCopyCommunication>
600 class InverseOperator final
601 : public Dune::Fem::Operator<typename LinearOperator::RangeFunctionType, typename LinearOperator::DomainFunctionType>
602 {
603 static_assert(std::is_same< typename LinearOperator::DomainFunctionType, typename LinearOperator::RangeFunctionType>::value,
604 "Domain function and range function must have the same type.");
605
606 public:
607 using LinearOperatorType = LinearOperator;
608 using OperatorType = LinearOperatorType;
609
610 using DiscreteFunctionType = typename LinearOperatorType::DomainFunctionType;
611 using DiscreteFunctionSpaceType = typename DiscreteFunctionType::DiscreteFunctionSpaceType;
612
613 using RealType = real_t<typename DiscreteFunctionType::RangeFieldType>;
614
615 using SolverParameterType = SolverParameter;
616
617 private:
618 using vector_type = VectorType<DiscreteFunctionType>;
619 using matrix_type = MatrixType<LinearOperatorType>;
620
621 public:
622 using CommunicationType = Communication<DiscreteFunctionSpaceType>;
623
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>;
627
628 using PreconditionerFactory = std::function<std::shared_ptr<PreconditionerType>(std::shared_ptr<AssembledLinearOperatorType>, const CommunicationType&, Symmetry, const SolverParameterType&)>;
629
630 InverseOperator (const SolverParameterType& parameter = {})
631 : InverseOperator(makePreconditioner<AssembledLinearOperatorType, CommunicationType>, parameter)
632 {}
633
634 InverseOperator (PreconditionerFactory factory, const SolverParameterType& parameter = {})
635 : preconditionerFactory_{std::move(factory)},
636 parameter_{std::make_shared<SolverParameterType>(parameter)}
637 {}
638
639 InverseOperator (const LinearOperatorType& op, const SolverParameterType& parameter = {})
640 : InverseOperator(parameter)
641 {
642 bind(op);
643 }
644
645 InverseOperator (const LinearOperatorType& op, PreconditionerFactory factory, const SolverParameterType& parameter = {})
646 : InverseOperator(std::move(factory), parameter)
647 {
648 bind(op);
649 }
650
651 void operator() (const DiscreteFunctionType& u, DiscreteFunctionType& w) const override
652 {
653 result_.clear();
654 vector_type b(u.dofVector().array());
655 if (parameter_->errorMeasure() == 0)
656 {
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_);
662 }
663 else
664 solver_->apply(w.dofVector().array(), b, result_);
665 }
666
667 void bind (const LinearOperatorType& op)
668 {
669 buildCommunication(op.domainSpace(), Dune::SolverCategory::overlapping, communication_);
670
671 linearOperator_.reset(new AssembledLinearOperatorType(op.exportMatrix(), *communication_));
672 scalarProduct_.reset(new Dune::OverlappingSchwarzScalarProduct<vector_type, CommunicationType>(*communication_));
673
674 // create preconditioner
675 preconditioner_ = preconditionerFactory_(linearOperator_, *communication_, symmetry, parameter());
676
677 // create linear solver
678 auto reduction = parameter().tolerance();
679 auto maxIterations = parameter().maxIterations();
680 auto verbosity = op.domainSpace().gridPart().comm().rank() == 0 ? parameter().verbosity() : 0;
681
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)))
690 {
691 case -2:
692 solver_.reset(new Dune::RestartedFlexibleGMResSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, parameter().restart(), maxIterations, verbosity));
693 break;
694
695 case -1:
696 solver_.reset(new Dune::RestartedFCGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity, parameter().fcgMmax()));
697 break;
698
699 case 0:
700 solver_.reset(new Dune::GeneralizedPCGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity, parameter().restart()));
701 break;
702
703 case 1:
704 solver_.reset(new Dune::CGSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
705 break;
706
707 case 2:
708 solver_.reset(new Dune::BiCGSTABSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
709 break;
710
711 case 3:
712 solver_.reset(new Dune::RestartedGMResSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, parameter().restart(), maxIterations, verbosity));
713 break;
714
715 case 4:
716 solver_.reset(new Dune::MINRESSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
717 break;
718
719 case 5:
720 solver_.reset(new Dune::GradientSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
721 break;
722
723 case 6:
724 solver_.reset(new Dune::LoopSolver<vector_type>(linearOperator_, scalarProduct_, preconditioner_, reduction, maxIterations, verbosity));
725 break;
726 }
727 }
728
729 void unbind ()
730 {
731 solver_.reset();
732 preconditioner_.reset();
733 scalarProduct_.reset();
734 linearOperator_.reset();
735 communication_.reset();
736 }
737
738 auto parameter () -> SolverParameterType& { return *parameter_; }
739 void setParamters (const SolverParameterType& parameter) { parameter_ = std::make_shared<SolverParameterType>(parameter); }
740
741 int iterations () const { return result_.iterations; }
742 bool converged () const { return result_.converged; }
743
744 void setMaxIterations (int maxIterations) { parameter().setMaxIterations(maxIterations); }
745
746 private:
747 PreconditionerFactory preconditionerFactory_;
748 std::shared_ptr<SolverParameterType> parameter_;
749
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_;
755
756 mutable Dune::InverseOperatorResult result_;
757 };
758
759
760
761 // OperatorBasedPreconditionerFactory
762 // ----------------------------------
763
764 template<class Operator, Symmetry symmetry = symmetric, template<class> class Communication = OwnerOverlapCopyCommunication>
765 class OperatorPreconditionerFactory
766 {
767 static_assert(std::is_same<typename Operator::DomainFunctionType, typename Operator::RangeFunctionType>::value,
768 "Domain function and range function must have the same type.");
769
770 using LinearOperatorType = typename Operator::JacobianOperatorType;
771 using DiscreteFunctionType = typename LinearOperatorType::DomainFunctionType;
772
773 using vector_type = VectorType<DiscreteFunctionType>;
774 using matrix_type = MatrixType<LinearOperatorType>;
775
776 public:
777 using DiscreteFunctionSpaceType = typename DiscreteFunctionType::DiscreteFunctionSpaceType;
778
779 using CommunicationType = Communication<DiscreteFunctionSpaceType>;
780
781 using PreconditionerType = Dune::Preconditioner<vector_type, vector_type>;
782
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)
788 {}
789
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>
794 {
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);
798 }
799
800 auto op () const -> const Operator& { return op_; }
801 auto op () -> Operator& { return op_; }
802
803 private:
804 Operator op_;
805 DiscreteFunctionType zero_;
806 mutable LinearOperatorType jacobian_;
807 };
808
809 } // namespace ISTL
810
811 } // namespace Fem
812
813} // namespace Dune
814
815#endif // #ifdef HAVE_DUNE_ISTL
816
817#endif // #ifndef DUNE_FEM_SOLVER_ISTL_HH
amg.hh
The AMG preconditioner.
bcrsmatrix.hh
Implementation of the BCRSMatrix class.
Dune::Amg::AMG
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:66
Dune::Amg::CoarsenCriterion
The criterion describing the stop criteria for the coarsening process.
Definition: matrixhierarchy.hh:283
Dune::Amg::Parameters
All parameters for AMG.
Definition: parameters.hh:416
Dune::Amg::SymmetricCriterion
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:519
Dune::Amg::UnSymmetricCriterion
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:539
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
Dune::BiCGSTABSolver
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:419
Dune::BlockPreconditioner
Block parallel preconditioner.
Definition: schwarz.hh:278
Dune::CGSolver
conjugate gradient method
Definition: solvers.hh:193
Dune::Fem::SparseRowMatrixObject::domainSpace
const DomainSpaceType & domainSpace() const
get domain space (i.e. space that builds the rows)
Definition: spmatrix.hh:633
Dune::Fem::TupleDiscreteFunction
forward declaration
Definition: discretefunction.hh:51
Dune::Fem::TupleDiscreteFunction::DiscreteFunctionSpaceType
TupleDiscreteFunctionSpace< typename DiscreteFunctions::DiscreteFunctionSpaceType ... > DiscreteFunctionSpaceType
type for the discrete function space this function lives in
Definition: discretefunction.hh:69
Dune::GeneralizedPCGSolver
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1307
Dune::GradientSolver
gradient method
Definition: solvers.hh:124
Dune::InvalidStateException
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Dune::LoopSolver
Preconditioned loop solver.
Definition: solvers.hh:59
Dune::MINRESSolver
Minimal Residual Method (MINRES)
Definition: solvers.hh:609
Dune::OverlappingSchwarzOperator
An overlapping Schwarz operator.
Definition: schwarz.hh:75
Dune::OverlappingSchwarzScalarProduct
Scalar product for overlapping Schwarz methods.
Definition: scalarproducts.hh:201
Dune::Preconditioner
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
Dune::RestartedFCGSolver
Accelerated flexible conjugate gradient method.
Definition: solvers.hh:1500
Dune::RestartedFlexibleGMResSolver
implements the Flexible Generalized Minimal Residual (FGMRes) method (right preconditioned)
Definition: solvers.hh:1139
Dune::RestartedGMResSolver
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:827
Dune::Richardson
Richardson preconditioner.
Definition: preconditioners.hh:878
Dune::ScalarProduct
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:52
Dune::SeqILDL
sequential ILDL preconditioner
Definition: preconditioners.hh:973
Dune::SeqILU
Sequential ILU preconditioner.
Definition: preconditioners.hh:697
Dune::SeqJac
The sequential jacobian preconditioner.
Definition: preconditioners.hh:413
Dune::SeqSOR
Sequential SOR preconditioner.
Definition: preconditioners.hh:262
Dune::SeqSSOR
Sequential SSOR preconditioner.
Definition: preconditioners.hh:142
exceptions.hh
A few common exception classes.
ftraits.hh
Type traits to determine the type of reals (when working with complex numbers)
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune::Hybrid::forEach
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
Dune::Amg::ConstructionTraits::Arguments
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
Dune::Amg::ConstructionTraits::construct
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:52
Dune::Amg::AccumulationMode
AccumulationMode
Identifiers for the different accumulation modes.
Definition: parameters.hh:231
Dune::Amg::DefaultSmootherArgs::relaxationFactor
RelaxationFactor relaxationFactor
The relaxation factor to use.
Definition: smoother.hh:51
Dune::Amg::DefaultSmootherArgs::iterations
int iterations
The number of iterations to perform.
Definition: smoother.hh:47
Dune
Dune namespace.
Definition: alignedallocator.hh:13
std
STL namespace.
operators.hh
Define general, extensible interface for operators. The available implementation wraps a matrix.
preconditioners.hh
Define general preconditioner interface.
solvers.hh
Implementations of the inverse operator interface.
Dune::Amg::DefaultSmootherArgs
The default class for the smoother arguments.
Definition: smoother.hh:38
Dune::Fem::Operator
abstract operator
Definition: operator.hh:34
Dune::Fem::Operator< DomainFunction, RangeFunction >::DomainFunctionType
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:36
Dune::Fem::SparseRowLinearOperator
SparseRowLinearOperator.
Definition: spoperator.hh:25
Dune::Fem::SparseRowLinearOperator::exportMatrix
MatrixType & exportMatrix() const
get reference to storage object
Definition: spmatrix.hh:655
Dune::InverseOperatorResult
Statistics about the application of an inverse operator.
Definition: solver.hh:50
Dune::MetaType
A type that refers to another type.
Definition: typelist.hh:33
Dune::SolverCategory::overlapping
@ overlapping
Category for overlapping solvers.
Definition: solvercategory.hh:29
typeutilities.hh
Utilities for type computations, constraining overloads, ...
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)