DUNE PDELab (2.7)

ovlpistlsolverbackend.hh
1// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=8 sw=2 sts=2:
3#ifndef DUNE_PDELAB_BACKEND_ISTL_OVLPISTLSOLVERBACKEND_HH
4#define DUNE_PDELAB_BACKEND_ISTL_OVLPISTLSOLVERBACKEND_HH
5
8
10#include <dune/istl/solvercategory.hh>
12#include <dune/istl/solvers.hh>
16#include <dune/istl/paamg/pinfo.hh>
17#include <dune/istl/io.hh>
18#include <dune/istl/superlu.hh>
19
20#include <dune/pdelab/constraints/common/constraints.hh>
21#include <dune/pdelab/gridfunctionspace/genericdatahandle.hh>
22#include <dune/pdelab/backend/istl/vector.hh>
23#include <dune/pdelab/backend/istl/bcrsmatrix.hh>
24#include <dune/pdelab/backend/istl/parallelhelper.hh>
25#include <dune/pdelab/backend/istl/seqistlsolverbackend.hh>
26
27namespace Dune {
28 namespace PDELab {
29
33
34 //========================================================
35 // Generic support for overlapping grids
36 // (need to be used with appropriate constraints)
37 //========================================================
38
39 // operator that resets result to zero at constrained DOFS
40 template<class CC, class M, class X, class Y>
41 class OverlappingOperator
42 : public Dune::AssembledLinearOperator<M,X,Y>
43 {
44 public:
46 typedef M matrix_type;
47 typedef X domain_type;
48 typedef Y range_type;
49 typedef typename X::ElementType field_type;
50
51 OverlappingOperator (const CC& cc_, const M& A)
52 : cc(cc_), _A_(A)
53 {}
54
56 virtual void apply (const domain_type& x, range_type& y) const override
57 {
58 using Backend::native;
59 native(_A_).mv(native(x),native(y));
61 }
62
64 virtual void applyscaleadd (field_type alpha, const domain_type& x, range_type& y) const override
65 {
66 using Backend::native;
67 native(_A_).usmv(alpha,native(x),native(y));
69 }
70
71 SolverCategory::Category category() const override
72 {
74 }
75
77 virtual const M& getmat () const override
78 {
79 return _A_;
80 }
81
82 private:
83 const CC& cc;
84 const M& _A_;
85 };
86
87 // new scalar product assuming at least overlap 1
88 // uses unique partitioning of nodes for parallelization
89 template<class GFS, class X>
90 class OverlappingScalarProduct
91 : public Dune::ScalarProduct<X>
92 {
93 public:
95 typedef X domain_type;
96 typedef typename X::ElementType field_type;
97
100 OverlappingScalarProduct (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_)
101 : gfs(gfs_), helper(helper_)
102 {}
103
104
109 virtual field_type dot (const X& x, const X& y) const override
110 {
111 // do local scalar product on unique partition
112 field_type sum = helper.disjointDot(x,y);
113
114 // do global communication
115 return gfs.gridView().comm().sum(sum);
116 }
117
121 virtual double norm (const X& x) const override
122 {
123 return sqrt(static_cast<double>(this->dot(x,x)));
124 }
125
126 SolverCategory::Category category() const override
127 {
129 }
130
131 private:
132 const GFS& gfs;
133 const ISTL::ParallelHelper<GFS>& helper;
134 };
135
136 // wrapped sequential preconditioner
137 template<class CC, class GFS, class P>
138 class OverlappingWrappedPreconditioner
139 : public Dune::Preconditioner<Dune::PDELab::Backend::Vector<GFS,typename P::domain_type::field_type>,
140 Dune::PDELab::Backend::Vector<GFS,typename P::range_type::field_type>>
141 {
142 public:
144 using domain_type = Dune::PDELab::Backend::Vector<GFS,typename P::domain_type::field_type>;
146 using range_type = Dune::PDELab::Backend::Vector<GFS,typename P::range_type::field_type>;
147
149 OverlappingWrappedPreconditioner (const GFS& gfs_, P& prec_, const CC& cc_,
150 const ISTL::ParallelHelper<GFS>& helper_)
151 : gfs(gfs_), prec(prec_), cc(cc_), helper(helper_)
152 {}
153
157 virtual void pre (domain_type& x, range_type& b) override
158 {
159 prec.pre(x,b);
160 }
161
165 virtual void apply (domain_type& v, const range_type& d) override
166 {
167 range_type dd(d);
168 set_constrained_dofs(cc,0.0,dd);
169 prec.apply(Backend::native(v),Backend::native(dd));
170 Dune::PDELab::AddDataHandle<GFS,domain_type> adddh(gfs,v);
171 if (gfs.gridView().comm().size()>1)
172 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
173 }
174
175 SolverCategory::Category category() const override
176 {
178 }
179
183 virtual void post (domain_type& x) override
184 {
185 prec.post(Backend::native(x));
186 }
187
188 private:
189 const GFS& gfs;
190 P& prec;
191 const CC& cc;
192 const ISTL::ParallelHelper<GFS>& helper;
193 };
194
195
196#if HAVE_SUITESPARSE_UMFPACK || DOXYGEN
197 // exact subdomain solves with UMFPack as preconditioner
198 template<class GFS, class M, class X, class Y>
199 class UMFPackSubdomainSolver : public Dune::Preconditioner<X,Y>
200 {
201 typedef Backend::Native<M> ISTLM;
202
203 public:
205 typedef X domain_type;
207 typedef Y range_type;
209 typedef typename X::ElementType field_type;
210
217 UMFPackSubdomainSolver (const GFS& gfs_, const M& A_)
218 : gfs(gfs_), solver(Backend::native(A_),false) // this does the decomposition
219 {}
220
224 virtual void pre (X& x, Y& b) override {}
225
229 virtual void apply (X& v, const Y& d) override
230 {
232 Y b(d); // need copy, since solver overwrites right hand side
233 solver.apply(Backend::native(v),Backend::native(b),stat);
234 if (gfs.gridView().comm().size()>1)
235 {
236 AddDataHandle<GFS,X> adddh(gfs,v);
237 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
238 }
239 }
240
241 SolverCategory::Category category() const override
242 {
244 }
245
249 virtual void post (X& x) override {}
250
251 private:
252 const GFS& gfs;
254 };
255#endif
256
257#if HAVE_SUPERLU
258 // exact subdomain solves with SuperLU as preconditioner
259 template<class GFS, class M, class X, class Y>
260 class SuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
261 {
262 typedef Backend::Native<M> ISTLM;
263
264 public:
266 typedef X domain_type;
268 typedef Y range_type;
270 typedef typename X::ElementType field_type;
271
278 SuperLUSubdomainSolver (const GFS& gfs_, const M& A_)
279 : gfs(gfs_), solver(Backend::native(A_),false) // this does the decomposition
280 {}
281
285 virtual void pre (X& x, Y& b) override {}
286
290 virtual void apply (X& v, const Y& d) override
291 {
293 Y b(d); // need copy, since solver overwrites right hand side
294 solver.apply(Backend::native(v),Backend::native(b),stat);
295 if (gfs.gridView().comm().size()>1)
296 {
297 AddDataHandle<GFS,X> adddh(gfs,v);
298 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
299 }
300 }
301
302 SolverCategory::Category category() const override
303 {
305 }
306
310 virtual void post (X& x) override {}
311
312 private:
313 const GFS& gfs;
315 };
316
317 // exact subdomain solves with SuperLU as preconditioner
318 template<class GFS, class M, class X, class Y>
319 class RestrictedSuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
320 {
321 typedef typename M::Container ISTLM;
322
323 public:
325 typedef X domain_type;
327 typedef Y range_type;
329 typedef typename X::ElementType field_type;
330
338 RestrictedSuperLUSubdomainSolver (const GFS& gfs_, const M& A_,
339 const ISTL::ParallelHelper<GFS>& helper_)
340 : gfs(gfs_), solver(Backend::native(A_),false), helper(helper_) // this does the decomposition
341 {}
342
346 virtual void pre (X& x, Y& b) override {}
347
351 virtual void apply (X& v, const Y& d) override
352 {
353 using Backend::native;
355 Y b(d); // need copy, since solver overwrites right hand side
356 solver.apply(native(v),native(b),stat);
357 if (gfs.gridView().comm().size()>1)
358 {
359 helper.maskForeignDOFs(native(v));
360 AddDataHandle<GFS,X> adddh(gfs,v);
361 gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
362 }
363 }
364
365 SolverCategory::Category category() const override
366 {
368 }
369
373 virtual void post (X& x) override {}
374
375 private:
376 const GFS& gfs;
378 const ISTL::ParallelHelper<GFS>& helper;
379 };
380#endif
381
382 template<typename GFS>
383 class OVLPScalarProductImplementation
384 {
385 public:
386 OVLPScalarProductImplementation(const GFS& gfs_)
387 : gfs(gfs_), helper(gfs_)
388 {}
389
394 template<typename X>
395 typename X::ElementType dot (const X& x, const X& y) const
396 {
397 // do local scalar product on unique partition
398 typename X::ElementType sum = helper.disjointDot(x,y);
399
400 // do global communication
401 return gfs.gridView().comm().sum(sum);
402 }
403
407 template<typename X>
408 typename Dune::template FieldTraits<typename X::ElementType >::real_type norm (const X& x) const
409 {
410 using namespace std;
411 return sqrt(static_cast<double>(this->dot(x,x)));
412 }
413
414 const ISTL::ParallelHelper<GFS>& parallelHelper() const
415 {
416 return helper;
417 }
418
419 // need also non-const version;
420 ISTL::ParallelHelper<GFS>& parallelHelper() // P.B.: needed for createIndexSetAndProjectForAMG
421 {
422 return helper;
423 }
424
425 private:
426 const GFS& gfs;
427 ISTL::ParallelHelper<GFS> helper;
428 };
429
430
431 template<typename GFS, typename X>
432 class OVLPScalarProduct
433 : public ScalarProduct<X>
434 {
435 public:
436 SolverCategory::Category category() const override
437 {
439 }
440
441 OVLPScalarProduct(const OVLPScalarProductImplementation<GFS>& implementation_)
442 : implementation(implementation_)
443 {}
444
445 virtual typename X::Container::field_type dot(const X& x, const X& y) const override
446 {
447 return implementation.dot(x,y);
448 }
449
450 virtual typename X::Container::field_type norm (const X& x) const override
451 {
452 using namespace std;
453 return sqrt(static_cast<double>(this->dot(x,x)));
454 }
455
456 private:
457 const OVLPScalarProductImplementation<GFS>& implementation;
458 };
459
460 template<class GFS, class C,
461 template<class,class,class,int> class Preconditioner,
462 template<class> class Solver>
463 class ISTLBackend_OVLP_Base
464 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
465 {
466 public:
475 ISTLBackend_OVLP_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
476 int steps_=5, int verbose_=1)
477 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), steps(steps_), verbose(verbose_)
478 {}
479
487 template<class M, class V, class W>
488 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
489 {
490 using Backend::Native;
491 using Backend::native;
492 typedef OverlappingOperator<C,M,V,W> POP;
493 POP pop(c,A);
494 typedef OVLPScalarProduct<GFS,V> PSP;
495 PSP psp(*this);
496 typedef Preconditioner<
497 Native<M>,
498 Native<V>,
499 Native<W>,
500 1
501 > SeqPrec;
502 SeqPrec seqprec(native(A),steps,1.0);
503 typedef OverlappingWrappedPreconditioner<C,GFS,SeqPrec> WPREC;
504 WPREC wprec(gfs,seqprec,c,this->parallelHelper());
505 int verb=0;
506 if (gfs.gridView().comm().rank()==0) verb=verbose;
507 Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
509 solver.apply(z,r,stat);
510 res.converged = stat.converged;
511 res.iterations = stat.iterations;
512 res.elapsed = stat.elapsed;
513 res.reduction = stat.reduction;
514 res.conv_rate = stat.conv_rate;
515 }
516 private:
517 const GFS& gfs;
518 const C& c;
519 unsigned maxiter;
520 int steps;
521 int verbose;
522 };
523
524 // Base class for ILU0 as preconditioner
525 template<class GFS, class C,
526 template<class> class Solver>
527 class ISTLBackend_OVLP_ILU0_Base
528 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
529 {
530 public:
538 ISTLBackend_OVLP_ILU0_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000, int verbose_=1)
539 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
540 {}
541
549 template<class M, class V, class W>
550 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
551 {
552 using Backend::Native;
553 using Backend::native;
554 typedef OverlappingOperator<C,M,V,W> POP;
555 POP pop(c,A);
556 typedef OVLPScalarProduct<GFS,V> PSP;
557 PSP psp(*this);
558 typedef SeqILU<
559 Native<M>,
560 Native<V>,
561 Native<W>,
562 1
563 > SeqPrec;
564 SeqPrec seqprec(native(A),1.0);
565 typedef OverlappingWrappedPreconditioner<C,GFS,SeqPrec> WPREC;
566 WPREC wprec(gfs,seqprec,c,this->parallelHelper());
567 int verb=0;
568 if (gfs.gridView().comm().rank()==0) verb=verbose;
569 Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
571 solver.apply(z,r,stat);
572 res.converged = stat.converged;
573 res.iterations = stat.iterations;
574 res.elapsed = stat.elapsed;
575 res.reduction = stat.reduction;
576 res.conv_rate = stat.conv_rate;
577 }
578 private:
579 const GFS& gfs;
580 const C& c;
581 unsigned maxiter;
582 int steps;
583 int verbose;
584 };
585
586 // Base class for ILUn as preconditioner
587 template<class GFS, class C,
588 template<class> class Solver>
589 class ISTLBackend_OVLP_ILUn_Base
590 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
591 {
592 public:
601 ISTLBackend_OVLP_ILUn_Base (const GFS& gfs_, const C& c_, int n_=1, unsigned maxiter_=5000, int verbose_=1)
602 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), n(n_), maxiter(maxiter_), verbose(verbose_)
603 {}
604
612 template<class M, class V, class W>
613 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
614 {
615 using Backend::Native;
616 using Backend::native;
617 typedef OverlappingOperator<C,M,V,W> POP;
618 POP pop(c,A);
619 typedef OVLPScalarProduct<GFS,V> PSP;
620 PSP psp(*this);
621 typedef SeqILU<
622 Native<M>,
623 Native<V>,
624 Native<W>,
625 1
626 > SeqPrec;
627 SeqPrec seqprec(native(A),n,1.0);
628 typedef OverlappingWrappedPreconditioner<C,GFS,SeqPrec> WPREC;
629 WPREC wprec(gfs,seqprec,c,this->parallelHelper());
630 int verb=0;
631 if (gfs.gridView().comm().rank()==0) verb=verbose;
632 Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
634 solver.apply(z,r,stat);
635 res.converged = stat.converged;
636 res.iterations = stat.iterations;
637 res.elapsed = stat.elapsed;
638 res.reduction = stat.reduction;
639 res.conv_rate = stat.conv_rate;
640 }
641 private:
642 const GFS& gfs;
643 const C& c;
644 int n;
645 unsigned maxiter;
646 int steps;
647 int verbose;
648 };
649
652
658 template<class GFS, class CC>
660 : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>
661 {
662 public:
671 ISTLBackend_OVLP_BCGS_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
672 int steps=5, int verbose=1)
673 : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>(gfs, cc, maxiter, steps, verbose)
674 {}
675 };
681 template<class GFS, class CC>
683 : public ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>
684 {
685 public:
693 ISTLBackend_OVLP_BCGS_ILU0 (const GFS& gfs, const CC& cc, unsigned maxiter=5000, int verbose=1)
694 : ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, maxiter, verbose)
695 {}
696 };
702 template<class GFS, class CC>
704 : public ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>
705 {
706 public:
715 ISTLBackend_OVLP_BCGS_ILUn (const GFS& gfs, const CC& cc, int n=1, unsigned maxiter=5000, int verbose=1)
716 : ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, n, maxiter, verbose)
717 {}
718 };
724 template<class GFS, class CC>
726 : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>
727 {
728 public:
737 ISTLBackend_OVLP_CG_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
738 int steps=5, int verbose=1)
739 : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>(gfs, cc, maxiter, steps, verbose)
740 {}
741 };
742
748 template<class GFS, class CC>
750 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
751 {
752 public:
760 ISTLBackend_OVLP_GMRES_ILU0 (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000, int verbose_=1,
761 int restart_ = 20)
762 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), cc(cc_), maxiter(maxiter_), verbose(verbose_),
763 restart(restart_)
764 {}
765
772 template<class M, class V, class W>
773 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
774 {
775 using Backend::Native;
776 using Backend::native;
777 typedef OverlappingOperator<CC,M,V,W> POP;
778 POP pop(cc,A);
779 typedef OVLPScalarProduct<GFS,V> PSP;
780 PSP psp(*this);
781 typedef SeqILU<
782 Native<M>,
783 Native<V>,
784 Native<W>,
785 1
786 > SeqPrec;
787 SeqPrec seqprec(native(A),1.0);
788 typedef OverlappingWrappedPreconditioner<CC,GFS,SeqPrec> WPREC;
789 WPREC wprec(gfs,seqprec,cc,this->parallelHelper());
790 int verb=0;
791 if (gfs.gridView().comm().rank()==0) verb=verbose;
792 RestartedGMResSolver<V> solver(pop,psp,wprec,reduction,restart,maxiter,verb);
794 solver.apply(z,r,stat);
795 res.converged = stat.converged;
796 res.iterations = stat.iterations;
797 res.elapsed = stat.elapsed;
798 res.reduction = stat.reduction;
799 res.conv_rate = stat.conv_rate;
800 }
801
802 private:
803 const GFS& gfs;
804 const CC& cc;
805 unsigned maxiter;
806 int steps;
807 int verbose;
808 int restart;
809 };
810
812
813 template<class GFS, class C, template<typename> class Solver>
814 class ISTLBackend_OVLP_SuperLU_Base
815 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
816 {
817 public:
825 ISTLBackend_OVLP_SuperLU_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
826 int verbose_=1)
827 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
828 {}
829
837 template<class M, class V, class W>
838 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
839 {
840 typedef OverlappingOperator<C,M,V,W> POP;
841 POP pop(c,A);
842 typedef OVLPScalarProduct<GFS,V> PSP;
843 PSP psp(*this);
844#if HAVE_SUPERLU
845 typedef SuperLUSubdomainSolver<GFS,M,V,W> PREC;
846 PREC prec(gfs,A);
847 int verb=0;
848 if (gfs.gridView().comm().rank()==0) verb=verbose;
849 Solver<V> solver(pop,psp,prec,reduction,maxiter,verb);
851 solver.apply(z,r,stat);
852 res.converged = stat.converged;
853 res.iterations = stat.iterations;
854 res.elapsed = stat.elapsed;
855 res.reduction = stat.reduction;
856 res.conv_rate = stat.conv_rate;
857#else
858 std::cout << "No superLU support, please install and configure it." << std::endl;
859#endif
860 }
861
862 private:
863 const GFS& gfs;
864 const C& c;
865 unsigned maxiter;
866 int verbose;
867 };
868
870
871 template<class GFS, class C, template<typename> class Solver>
872 class ISTLBackend_OVLP_UMFPack_Base
873 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
874 {
875 public:
883 ISTLBackend_OVLP_UMFPack_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
884 int verbose_=1)
885 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
886 {}
887
895 template<class M, class V, class W>
896 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
897 {
898 typedef OverlappingOperator<C,M,V,W> POP;
899 POP pop(c,A);
900 typedef OVLPScalarProduct<GFS,V> PSP;
901 PSP psp(*this);
902#if HAVE_SUITESPARSE_UMFPACK || DOXYGEN
903 typedef UMFPackSubdomainSolver<GFS,M,V,W> PREC;
904 PREC prec(gfs,A);
905 int verb=0;
906 if (gfs.gridView().comm().rank()==0) verb=verbose;
907 Solver<V> solver(pop,psp,prec,reduction,maxiter,verb);
909 solver.apply(z,r,stat);
910 res.converged = stat.converged;
911 res.iterations = stat.iterations;
912 res.elapsed = stat.elapsed;
913 res.reduction = stat.reduction;
914 res.conv_rate = stat.conv_rate;
915#else
916 std::cout << "No UMFPack support, please install and configure it." << std::endl;
917#endif
918 }
919
920 private:
921 const GFS& gfs;
922 const C& c;
923 unsigned maxiter;
924 int verbose;
925 };
926
929
934 template<class GFS, class CC>
936 : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>
937 {
938 public:
939
947 ISTLBackend_OVLP_BCGS_SuperLU (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000,
948 int verbose_=1)
949 : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs_,cc_,maxiter_,verbose_)
950 {}
951 };
952
958 template<class GFS, class CC>
960 : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>
961 {
962 public:
963
971 ISTLBackend_OVLP_CG_SuperLU (const GFS& gfs_, const CC& cc_,
972 unsigned maxiter_=5000,
973 int verbose_=1)
974 : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
975 {}
976 };
977
983 template<class GFS, class CC>
985 : public ISTLBackend_OVLP_UMFPack_Base<GFS,CC,Dune::CGSolver>
986 {
987 public:
988
996 ISTLBackend_OVLP_CG_UMFPack (const GFS& gfs_, const CC& cc_,
997 unsigned maxiter_=5000,
998 int verbose_=1)
999 : ISTLBackend_OVLP_UMFPack_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
1000 {}
1001 };
1002
1003
1007 template<class GFS>
1009 : public LinearResultStorage
1010 {
1011 public:
1016 explicit ISTLBackend_OVLP_ExplicitDiagonal (const GFS& gfs_)
1017 : gfs(gfs_)
1018 {}
1019
1021 : gfs(other_.gfs)
1022 {}
1023
1028 template<class V>
1029 typename V::ElementType norm(const V& v) const
1030 {
1031 static_assert
1033 "ISTLBackend_OVLP_ExplicitDiagonal::norm() should not be "
1034 "neccessary, so we skipped the implementation. If you have a "
1035 "scenario where you need it, please implement it or report back to "
1036 "us.");
1037 }
1038
1046 template<class M, class V, class W>
1047 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
1048 {
1049 using Backend::Native;
1050 using Backend::native;
1052 Native<M>,
1053 Native<V>,
1054 Native<W>
1055 > jac(native(A),1,1.0);
1056 jac.pre(native(z),native(r));
1057 jac.apply(native(z),native(r));
1058 jac.post(native(z));
1059 if (gfs.gridView().comm().size()>1)
1060 {
1061 CopyDataHandle<GFS,V> copydh(gfs,z);
1062 gfs.gridView().communicate(copydh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
1063 }
1064 res.converged = true;
1065 res.iterations = 1;
1066 res.elapsed = 0.0;
1067 res.reduction = static_cast<double>(reduction);
1068 res.conv_rate = static_cast<double>(reduction); // pow(reduction,1.0/1)
1069 }
1070
1071 private:
1072 const GFS& gfs;
1073 };
1075
1076 template<class GO, int s, template<class,class,class,int> class Preconditioner,
1077 template<class> class Solver>
1078 class ISTLBackend_AMG : public LinearResultStorage
1079 {
1080 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1081 typedef ISTL::ParallelHelper<GFS> PHELPER;
1082 typedef typename GO::Traits::Jacobian M;
1083 typedef Backend::Native<M> MatrixType;
1084 typedef typename GO::Traits::Domain V;
1085 typedef Backend::Native<V> VectorType;
1086 typedef typename ISTL::CommSelector<s,Dune::MPIHelper::isFake>::type Comm;
1087#if HAVE_MPI
1091#else
1094#endif
1095 typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
1097
1098 typedef typename V::ElementType RF;
1099
1100 public:
1101
1106
1107 public:
1108 ISTLBackend_AMG(const GFS& gfs_, unsigned maxiter_=5000,
1109 int verbose_=1, bool reuse_=false,
1110 bool usesuperlu_=true)
1111 : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), params(15,2000),
1112 verbose(verbose_), reuse(reuse_), firstapply(true),
1113 usesuperlu(usesuperlu_)
1114 {
1115 params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
1116 params.setDebugLevel(verbose_);
1117#if !HAVE_SUPERLU
1118 if (gfs.gridView().comm().rank() == 0 && usesuperlu == true)
1119 {
1120 std::cout << "WARNING: You are using AMG without SuperLU!"
1121 << " Please consider installing SuperLU,"
1122 << " or set the usesuperlu flag to false"
1123 << " to suppress this warning." << std::endl;
1124 }
1125#endif
1126 }
1127
1132 void setParameters(const Parameters& params_)
1133 {
1134 params = params_;
1135 }
1136
1144 const Parameters& parameters() const
1145 {
1146 return params;
1147 }
1148
1150 void setReuse(bool reuse_)
1151 {
1152 reuse = reuse_;
1153 }
1154
1156 bool getReuse() const
1157 {
1158 return reuse;
1159 }
1160
1165 typename V::ElementType norm (const V& v) const
1166 {
1167 typedef OverlappingScalarProduct<GFS,V> PSP;
1168 PSP psp(gfs,phelper);
1169 return psp.norm(v);
1170 }
1171
1179 void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
1180 {
1181 Timer watch;
1182 Comm oocc(gfs.gridView().comm());
1183 MatrixType& mat=Backend::native(A);
1185 Dune::Amg::FirstDiagonal> > Criterion;
1186#if HAVE_MPI
1187 phelper.createIndexSetAndProjectForAMG(A, oocc);
1188 Operator oop(mat, oocc);
1190#else
1191 Operator oop(mat);
1193#endif
1194 SmootherArgs smootherArgs;
1195 smootherArgs.iterations = 1;
1196 smootherArgs.relaxationFactor = 1;
1197 Criterion criterion(params);
1198 stats.tprepare=watch.elapsed();
1199 watch.reset();
1200
1201 int verb=0;
1202 if (gfs.gridView().comm().rank()==0) verb=verbose;
1203 //only construct a new AMG if the matrix changes
1204 if (reuse==false || firstapply==true){
1205 amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
1206 firstapply = false;
1207 stats.tsetup = watch.elapsed();
1208 stats.levels = amg->maxlevels();
1209 stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
1210 }
1211 watch.reset();
1212 Solver<VectorType> solver(oop,sp,*amg,RF(reduction),maxiter,verb);
1214
1215 solver.apply(Backend::native(z),Backend::native(r),stat);
1216 stats.tsolve= watch.elapsed();
1217 res.converged = stat.converged;
1218 res.iterations = stat.iterations;
1219 res.elapsed = stat.elapsed;
1220 res.reduction = stat.reduction;
1221 res.conv_rate = stat.conv_rate;
1222 }
1223
1228 const ISTLAMGStatistics& statistics() const
1229 {
1230 return stats;
1231 }
1232
1233 private:
1234 const GFS& gfs;
1235 PHELPER phelper;
1236 unsigned maxiter;
1237 Parameters params;
1238 int verbose;
1239 bool reuse;
1240 bool firstapply;
1241 bool usesuperlu;
1242 std::shared_ptr<AMG> amg;
1243 ISTLAMGStatistics stats;
1244 };
1245
1248
1255 template<class GO, int s=96>
1257 : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1258 {
1259 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1260 public:
1270 ISTLBackend_CG_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1271 int verbose_=1, bool reuse_=false,
1272 bool usesuperlu_=true)
1273 : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1274 (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1275 {}
1276 };
1277
1284 template<class GO, int s=96>
1286 : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1287 {
1288 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1289 public:
1299 ISTLBackend_BCGS_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1300 int verbose_=1, bool reuse_=false,
1301 bool usesuperlu_=true)
1302 : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1303 (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1304 {}
1305 };
1306
1313 template<class GO, int s=96>
1315 : public ISTLBackend_AMG<GO, s, Dune::SeqILU, Dune::BiCGSTABSolver>
1316 {
1317 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1318 public:
1328 ISTLBackend_BCGS_AMG_ILU0(const GFS& gfs_, unsigned maxiter_=5000,
1329 int verbose_=1, bool reuse_=false,
1330 bool usesuperlu_=true)
1331 : ISTLBackend_AMG<GO, s, Dune::SeqILU, Dune::BiCGSTABSolver>
1332 (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1333 {}
1334 };
1335
1337
1339
1340 } // namespace PDELab
1341} // namespace Dune
1342
1343#endif // DUNE_PDELAB_BACKEND_ISTL_OVLPISTLSOLVERBACKEND_HH
The AMG preconditioner.
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:62
The criterion describing the stop criteria for the coarsening process.
Definition: matrixhierarchy.hh:283
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:451
All parameters for AMG.
Definition: parameters.hh:391
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:515
A linear operator exporting itself in matrix form.
Definition: operators.hh:107
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:417
Block parallel preconditioner.
Definition: schwarz.hh:273
conjugate gradient method
Definition: solvers.hh:189
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:135
An overlapping Schwarz operator.
Definition: schwarz.hh:76
Scalar product for overlapping Schwarz methods.
Definition: scalarproducts.hh:183
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by ILU0.
Definition: ovlpistlsolverbackend.hh:1316
ISTLBackend_BCGS_AMG_ILU0(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1328
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR.
Definition: ovlpistlsolverbackend.hh:1287
ISTLBackend_BCGS_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1299
Overlapping parallel conjugate gradient solver preconditioned with AMG smoothed by SSOR.
Definition: ovlpistlsolverbackend.hh:1258
ISTLBackend_CG_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1270
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:684
ISTLBackend_OVLP_BCGS_ILU0(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:693
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:705
ISTLBackend_OVLP_BCGS_ILUn(const GFS &gfs, const CC &cc, int n=1, unsigned maxiter=5000, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:715
Overlapping parallel BiCGStab solver with SSOR preconditioner.
Definition: ovlpistlsolverbackend.hh:661
ISTLBackend_OVLP_BCGS_SSORk(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int steps=5, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:671
Overlapping parallel BiCGStab solver with SuperLU preconditioner.
Definition: ovlpistlsolverbackend.hh:937
ISTLBackend_OVLP_BCGS_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:947
Overlapping parallel CGS solver with SSOR preconditioner.
Definition: ovlpistlsolverbackend.hh:727
ISTLBackend_OVLP_CG_SSORk(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int steps=5, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:737
Overlapping parallel CG solver with SuperLU preconditioner.
Definition: ovlpistlsolverbackend.hh:961
ISTLBackend_OVLP_CG_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:971
Overlapping parallel CG solver with UMFPack preconditioner.
Definition: ovlpistlsolverbackend.hh:986
ISTLBackend_OVLP_CG_UMFPack(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:996
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: ovlpistlsolverbackend.hh:1010
ISTLBackend_OVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:1016
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: ovlpistlsolverbackend.hh:1029
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename W::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:1047
Overlapping parallel restarted GMRes solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:751
ISTLBackend_OVLP_GMRES_ILU0(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1, int restart_=20)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:760
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:773
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
virtual void pre(Dune::PDELab::Backend::Vector< GFS, P::domain_type::field_type > &x, Dune::PDELab::Backend::Vector< GFS, P::range_type::field_type > &b)=0
Prepare the preconditioner.
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:812
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:895
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:48
Sequential ILU preconditioner.
Definition: preconditioners.hh:500
The sequential jacobian preconditioner.
Definition: preconditioners.hh:390
Sequential SSOR preconditioner.
Definition: preconditioners.hh:138
Default implementation for the scalar case.
Definition: scalarproducts.hh:153
Definition: recipe-operator-splitting.cc:108
Definition of the DUNE_DEPRECATED macro for the case that config.h is not available.
Some generic functions for pretty printing vectors and matrices.
Implementations of the inverse operator interface.
auto dot(const A &a, const B &b) -> typename std::enable_if<!IsVector< A >::value &&!std::is_same< typename FieldTraits< A >::field_type, typename FieldTraits< A >::real_type > ::value, decltype(conj(a) *b)>::type
computes the dot product for fundamental data types according to Petsc's VectDot function: dot(a,...
Definition: dotproduct.hh:40
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:46
@ ForwardCommunication
communicate as given in InterfaceType
Definition: gridenums.hh:169
@ InteriorBorder_All_Interface
send interior and border, receive all entities
Definition: gridenums.hh:86
@ All_All_Interface
send all and receive all entities
Definition: gridenums.hh:89
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
Helpers for dealing with MPI.
Dune namespace.
Definition: alignedallocator.hh:14
STL namespace.
Define general, extensible interface for operators. The available implementation wraps a matrix.
Classes providing communication interfaces for overlapping Schwarz methods.
Define general preconditioner interface.
Define base class for scalar product and norm.
template which always yields a false value
Definition: typetraits.hh:126
The default class for the smoother arguments.
Definition: smoother.hh:37
Statistics about the application of an inverse operator.
Definition: solver.hh:46
double elapsed
Elapsed time in seconds.
Definition: solver.hh:80
int iterations
Number of iterations.
Definition: solver.hh:65
double reduction
Reduction achieved: .
Definition: solver.hh:68
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:74
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Category
Definition: solvercategory.hh:21
@ overlapping
Category for overlapping solvers.
Definition: solvercategory.hh:27
Classes for using SuperLU with ISTL matrices.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)