DUNE PDELab (git)

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 using std::sqrt;
124 return sqrt(static_cast<double>(this->dot(x,x)));
125 }
126
127 SolverCategory::Category category() const override
128 {
130 }
131
132 private:
133 const GFS& gfs;
134 const ISTL::ParallelHelper<GFS>& helper;
135 };
136
137 // wrapped sequential preconditioner
138 template<class CC, class GFS, class P>
139 class OverlappingWrappedPreconditioner
140 : public Dune::Preconditioner<Dune::PDELab::Backend::Vector<GFS,typename P::domain_type::field_type>,
141 Dune::PDELab::Backend::Vector<GFS,typename P::range_type::field_type>>
142 {
143 public:
145 using domain_type = Dune::PDELab::Backend::Vector<GFS,typename P::domain_type::field_type>;
147 using range_type = Dune::PDELab::Backend::Vector<GFS,typename P::range_type::field_type>;
148
150 OverlappingWrappedPreconditioner (const GFS& gfs_, P& prec_, const CC& cc_,
151 const ISTL::ParallelHelper<GFS>& helper_)
152 : gfs(gfs_), prec(prec_), cc(cc_), helper(helper_)
153 {}
154
158 virtual void pre (domain_type& x, range_type& b) override
159 {
160 prec.pre(x,b);
161 }
162
166 virtual void apply (domain_type& v, const range_type& d) override
167 {
168 range_type dd(d);
169 set_constrained_dofs(cc,0.0,dd);
170 prec.apply(Backend::native(v),Backend::native(dd));
171 Dune::PDELab::AddDataHandle<GFS,domain_type> adddh(gfs,v);
172 if (gfs.gridView().comm().size()>1)
173 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
174 }
175
176 SolverCategory::Category category() const override
177 {
179 }
180
184 virtual void post (domain_type& x) override
185 {
186 prec.post(Backend::native(x));
187 }
188
189 private:
190 const GFS& gfs;
191 P& prec;
192 const CC& cc;
193 const ISTL::ParallelHelper<GFS>& helper;
194 };
195
196
197#if HAVE_SUITESPARSE_UMFPACK || DOXYGEN
198 // exact subdomain solves with UMFPack as preconditioner
199 template<class GFS, class M, class X, class Y>
200 class UMFPackSubdomainSolver : public Dune::Preconditioner<X,Y>
201 {
202 typedef Backend::Native<M> ISTLM;
203
204 public:
206 typedef X domain_type;
208 typedef Y range_type;
210 typedef typename X::ElementType field_type;
211
218 UMFPackSubdomainSolver (const GFS& gfs_, const M& A_)
219 : gfs(gfs_), solver(Backend::native(A_),false) // this does the decomposition
220 {}
221
225 virtual void pre (X& x, Y& b) override {}
226
230 virtual void apply (X& v, const Y& d) override
231 {
233 Y b(d); // need copy, since solver overwrites right hand side
234 solver.apply(Backend::native(v),Backend::native(b),stat);
235 if (gfs.gridView().comm().size()>1)
236 {
237 AddDataHandle<GFS,X> adddh(gfs,v);
238 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
239 }
240 }
241
242 SolverCategory::Category category() const override
243 {
245 }
246
250 virtual void post (X& x) override {}
251
252 private:
253 const GFS& gfs;
255 };
256#endif
257
258#if HAVE_SUPERLU
259 // exact subdomain solves with SuperLU as preconditioner
260 template<class GFS, class M, class X, class Y>
261 class SuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
262 {
263 typedef Backend::Native<M> ISTLM;
264
265 public:
267 typedef X domain_type;
269 typedef Y range_type;
271 typedef typename X::ElementType field_type;
272
279 SuperLUSubdomainSolver (const GFS& gfs_, const M& A_)
280 : gfs(gfs_), solver(Backend::native(A_),false) // this does the decomposition
281 {}
282
286 virtual void pre (X& x, Y& b) override {}
287
291 virtual void apply (X& v, const Y& d) override
292 {
294 Y b(d); // need copy, since solver overwrites right hand side
295 solver.apply(Backend::native(v),Backend::native(b),stat);
296 if (gfs.gridView().comm().size()>1)
297 {
298 AddDataHandle<GFS,X> adddh(gfs,v);
299 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
300 }
301 }
302
303 SolverCategory::Category category() const override
304 {
306 }
307
311 virtual void post (X& x) override {}
312
313 private:
314 const GFS& gfs;
316 };
317
318 // exact subdomain solves with SuperLU as preconditioner
319 template<class GFS, class M, class X, class Y>
320 class RestrictedSuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
321 {
322 typedef typename M::Container ISTLM;
323
324 public:
326 typedef X domain_type;
328 typedef Y range_type;
330 typedef typename X::ElementType field_type;
331
339 RestrictedSuperLUSubdomainSolver (const GFS& gfs_, const M& A_,
340 const ISTL::ParallelHelper<GFS>& helper_)
341 : gfs(gfs_), solver(Backend::native(A_),false), helper(helper_) // this does the decomposition
342 {}
343
347 virtual void pre (X& x, Y& b) override {}
348
352 virtual void apply (X& v, const Y& d) override
353 {
354 using Backend::native;
356 Y b(d); // need copy, since solver overwrites right hand side
357 solver.apply(native(v),native(b),stat);
358 if (gfs.gridView().comm().size()>1)
359 {
360 helper.maskForeignDOFs(native(v));
361 AddDataHandle<GFS,X> adddh(gfs,v);
362 gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
363 }
364 }
365
366 SolverCategory::Category category() const override
367 {
369 }
370
374 virtual void post (X& x) override {}
375
376 private:
377 const GFS& gfs;
379 const ISTL::ParallelHelper<GFS>& helper;
380 };
381#endif
382
383 template<typename GFS>
384 class OVLPScalarProductImplementation
385 {
386 public:
387 OVLPScalarProductImplementation(const GFS& gfs_)
388 : gfs(gfs_), helper(gfs_)
389 {}
390
395 template<typename X>
396 typename X::ElementType dot (const X& x, const X& y) const
397 {
398 // do local scalar product on unique partition
399 typename X::ElementType sum = helper.disjointDot(x,y);
400
401 // do global communication
402 return gfs.gridView().comm().sum(sum);
403 }
404
408 template<typename X>
409 typename Dune::template FieldTraits<typename X::ElementType >::real_type norm (const X& x) const
410 {
411 using std::sqrt;
412 return sqrt(static_cast<double>(this->dot(x,x)));
413 }
414
415 const ISTL::ParallelHelper<GFS>& parallelHelper() const
416 {
417 return helper;
418 }
419
420 // need also non-const version;
421 ISTL::ParallelHelper<GFS>& parallelHelper() // P.B.: needed for createIndexSetAndProjectForAMG
422 {
423 return helper;
424 }
425
426 private:
427 const GFS& gfs;
428 ISTL::ParallelHelper<GFS> helper;
429 };
430
431
432 template<typename GFS, typename X>
433 class OVLPScalarProduct
434 : public ScalarProduct<X>
435 {
436 public:
437 SolverCategory::Category category() const override
438 {
440 }
441
442 OVLPScalarProduct(const OVLPScalarProductImplementation<GFS>& implementation_)
443 : implementation(implementation_)
444 {}
445
446 virtual typename X::Container::field_type dot(const X& x, const X& y) const override
447 {
448 return implementation.dot(x,y);
449 }
450
451 virtual typename X::Container::field_type norm (const X& x) const override
452 {
453 using std::sqrt;
454 return sqrt(static_cast<double>(this->dot(x,x)));
455 }
456
457 private:
458 const OVLPScalarProductImplementation<GFS>& implementation;
459 };
460
461 template<class GFS, class C,
462 template<class,class,class,int> class Preconditioner,
463 template<class> class Solver>
464 class ISTLBackend_OVLP_Base
465 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
466 {
467 public:
476 ISTLBackend_OVLP_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
477 int steps_=5, int verbose_=1)
478 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), steps(steps_), verbose(verbose_)
479 {}
480
488 template<class M, class V, class W>
489 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
490 {
491 using Backend::Native;
492 using Backend::native;
493 typedef OverlappingOperator<C,M,V,W> POP;
494 POP pop(c,A);
495 typedef OVLPScalarProduct<GFS,V> PSP;
496 PSP psp(*this);
497 typedef Preconditioner<
498 Native<M>,
499 Native<V>,
500 Native<W>,
501 1
502 > SeqPrec;
503 SeqPrec seqprec(native(A),steps,1.0);
504 typedef OverlappingWrappedPreconditioner<C,GFS,SeqPrec> WPREC;
505 WPREC wprec(gfs,seqprec,c,this->parallelHelper());
506 int verb=0;
507 if (gfs.gridView().comm().rank()==0) verb=verbose;
508 Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
510 solver.apply(z,r,stat);
511 res.converged = stat.converged;
512 res.iterations = stat.iterations;
513 res.elapsed = stat.elapsed;
514 res.reduction = stat.reduction;
515 res.conv_rate = stat.conv_rate;
516 }
517 private:
518 const GFS& gfs;
519 const C& c;
520 unsigned maxiter;
521 int steps;
522 int verbose;
523 };
524
525 // Base class for ILU0 as preconditioner
526 template<class GFS, class C,
527 template<class> class Solver>
528 class ISTLBackend_OVLP_ILU0_Base
529 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
530 {
531 public:
539 ISTLBackend_OVLP_ILU0_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000, int verbose_=1)
540 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
541 {}
542
550 template<class M, class V, class W>
551 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
552 {
553 using Backend::Native;
554 using Backend::native;
555 typedef OverlappingOperator<C,M,V,W> POP;
556 POP pop(c,A);
557 typedef OVLPScalarProduct<GFS,V> PSP;
558 PSP psp(*this);
559 typedef SeqILU<
560 Native<M>,
561 Native<V>,
562 Native<W>,
563 1
564 > SeqPrec;
565 SeqPrec seqprec(native(A),1.0);
566 typedef OverlappingWrappedPreconditioner<C,GFS,SeqPrec> WPREC;
567 WPREC wprec(gfs,seqprec,c,this->parallelHelper());
568 int verb=0;
569 if (gfs.gridView().comm().rank()==0) verb=verbose;
570 Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
572 solver.apply(z,r,stat);
573 res.converged = stat.converged;
574 res.iterations = stat.iterations;
575 res.elapsed = stat.elapsed;
576 res.reduction = stat.reduction;
577 res.conv_rate = stat.conv_rate;
578 }
579 private:
580 const GFS& gfs;
581 const C& c;
582 unsigned maxiter;
583 int steps;
584 int verbose;
585 };
586
587 // Base class for ILUn as preconditioner
588 template<class GFS, class C,
589 template<class> class Solver>
590 class ISTLBackend_OVLP_ILUn_Base
591 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
592 {
593 public:
602 ISTLBackend_OVLP_ILUn_Base (const GFS& gfs_, const C& c_, int n_=1, unsigned maxiter_=5000, int verbose_=1)
603 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), n(n_), maxiter(maxiter_), verbose(verbose_)
604 {}
605
613 template<class M, class V, class W>
614 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
615 {
616 using Backend::Native;
617 using Backend::native;
618 typedef OverlappingOperator<C,M,V,W> POP;
619 POP pop(c,A);
620 typedef OVLPScalarProduct<GFS,V> PSP;
621 PSP psp(*this);
622 typedef SeqILU<
623 Native<M>,
624 Native<V>,
625 Native<W>,
626 1
627 > SeqPrec;
628 SeqPrec seqprec(native(A),n,1.0);
629 typedef OverlappingWrappedPreconditioner<C,GFS,SeqPrec> WPREC;
630 WPREC wprec(gfs,seqprec,c,this->parallelHelper());
631 int verb=0;
632 if (gfs.gridView().comm().rank()==0) verb=verbose;
633 Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
635 solver.apply(z,r,stat);
636 res.converged = stat.converged;
637 res.iterations = stat.iterations;
638 res.elapsed = stat.elapsed;
639 res.reduction = stat.reduction;
640 res.conv_rate = stat.conv_rate;
641 }
642 private:
643 const GFS& gfs;
644 const C& c;
645 int n;
646 unsigned maxiter;
647 int steps;
648 int verbose;
649 };
650
653
659 template<class GFS, class CC>
661 : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>
662 {
663 public:
672 ISTLBackend_OVLP_BCGS_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
673 int steps=5, int verbose=1)
674 : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>(gfs, cc, maxiter, steps, verbose)
675 {}
676 };
682 template<class GFS, class CC>
684 : public ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>
685 {
686 public:
694 ISTLBackend_OVLP_BCGS_ILU0 (const GFS& gfs, const CC& cc, unsigned maxiter=5000, int verbose=1)
695 : ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, maxiter, verbose)
696 {}
697 };
703 template<class GFS, class CC>
705 : public ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>
706 {
707 public:
716 ISTLBackend_OVLP_BCGS_ILUn (const GFS& gfs, const CC& cc, int n=1, unsigned maxiter=5000, int verbose=1)
717 : ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, n, maxiter, verbose)
718 {}
719 };
725 template<class GFS, class CC>
727 : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>
728 {
729 public:
738 ISTLBackend_OVLP_CG_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
739 int steps=5, int verbose=1)
740 : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>(gfs, cc, maxiter, steps, verbose)
741 {}
742 };
743
749 template<class GFS, class CC>
751 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
752 {
753 public:
761 ISTLBackend_OVLP_GMRES_ILU0 (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000, int verbose_=1,
762 int restart_ = 20)
763 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), cc(cc_), maxiter(maxiter_), verbose(verbose_),
764 restart(restart_)
765 {}
766
773 template<class M, class V, class W>
774 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
775 {
776 using Backend::Native;
777 using Backend::native;
778 typedef OverlappingOperator<CC,M,V,W> POP;
779 POP pop(cc,A);
780 typedef OVLPScalarProduct<GFS,V> PSP;
781 PSP psp(*this);
782 typedef SeqILU<
783 Native<M>,
784 Native<V>,
785 Native<W>,
786 1
787 > SeqPrec;
788 SeqPrec seqprec(native(A),1.0);
789 typedef OverlappingWrappedPreconditioner<CC,GFS,SeqPrec> WPREC;
790 WPREC wprec(gfs,seqprec,cc,this->parallelHelper());
791 int verb=0;
792 if (gfs.gridView().comm().rank()==0) verb=verbose;
793 RestartedGMResSolver<V> solver(pop,psp,wprec,reduction,restart,maxiter,verb);
795 solver.apply(z,r,stat);
796 res.converged = stat.converged;
797 res.iterations = stat.iterations;
798 res.elapsed = stat.elapsed;
799 res.reduction = stat.reduction;
800 res.conv_rate = stat.conv_rate;
801 }
802
803 private:
804 const GFS& gfs;
805 const CC& cc;
806 unsigned maxiter;
807 int steps;
808 int verbose;
809 int restart;
810 };
811
813
814 template<class GFS, class C, template<typename> class Solver>
815 class ISTLBackend_OVLP_SuperLU_Base
816 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
817 {
818 public:
826 ISTLBackend_OVLP_SuperLU_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
827 int verbose_=1)
828 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
829 {}
830
838 template<class M, class V, class W>
839 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
840 {
841 typedef OverlappingOperator<C,M,V,W> POP;
842 POP pop(c,A);
843 typedef OVLPScalarProduct<GFS,V> PSP;
844 PSP psp(*this);
845#if HAVE_SUPERLU
846 typedef SuperLUSubdomainSolver<GFS,M,V,W> PREC;
847 PREC prec(gfs,A);
848 int verb=0;
849 if (gfs.gridView().comm().rank()==0) verb=verbose;
850 Solver<V> solver(pop,psp,prec,reduction,maxiter,verb);
852 solver.apply(z,r,stat);
853 res.converged = stat.converged;
854 res.iterations = stat.iterations;
855 res.elapsed = stat.elapsed;
856 res.reduction = stat.reduction;
857 res.conv_rate = stat.conv_rate;
858#else
859 std::cout << "No superLU support, please install and configure it." << std::endl;
860#endif
861 }
862
863 private:
864 const GFS& gfs;
865 const C& c;
866 unsigned maxiter;
867 int verbose;
868 };
869
871
872 template<class GFS, class C, template<typename> class Solver>
873 class ISTLBackend_OVLP_UMFPack_Base
874 : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
875 {
876 public:
884 ISTLBackend_OVLP_UMFPack_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
885 int verbose_=1)
886 : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
887 {}
888
896 template<class M, class V, class W>
897 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
898 {
899 typedef OverlappingOperator<C,M,V,W> POP;
900 POP pop(c,A);
901 typedef OVLPScalarProduct<GFS,V> PSP;
902 PSP psp(*this);
903#if HAVE_SUITESPARSE_UMFPACK || DOXYGEN
904 typedef UMFPackSubdomainSolver<GFS,M,V,W> PREC;
905 PREC prec(gfs,A);
906 int verb=0;
907 if (gfs.gridView().comm().rank()==0) verb=verbose;
908 Solver<V> solver(pop,psp,prec,reduction,maxiter,verb);
910 solver.apply(z,r,stat);
911 res.converged = stat.converged;
912 res.iterations = stat.iterations;
913 res.elapsed = stat.elapsed;
914 res.reduction = stat.reduction;
915 res.conv_rate = stat.conv_rate;
916#else
917 std::cout << "No UMFPack support, please install and configure it." << std::endl;
918#endif
919 }
920
921 private:
922 const GFS& gfs;
923 const C& c;
924 unsigned maxiter;
925 int verbose;
926 };
927
930
935 template<class GFS, class CC>
937 : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>
938 {
939 public:
940
948 ISTLBackend_OVLP_BCGS_SuperLU (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000,
949 int verbose_=1)
950 : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs_,cc_,maxiter_,verbose_)
951 {}
952 };
953
959 template<class GFS, class CC>
961 : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>
962 {
963 public:
964
972 ISTLBackend_OVLP_CG_SuperLU (const GFS& gfs_, const CC& cc_,
973 unsigned maxiter_=5000,
974 int verbose_=1)
975 : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
976 {}
977 };
978
984 template<class GFS, class CC>
986 : public ISTLBackend_OVLP_UMFPack_Base<GFS,CC,Dune::CGSolver>
987 {
988 public:
989
997 ISTLBackend_OVLP_CG_UMFPack (const GFS& gfs_, const CC& cc_,
998 unsigned maxiter_=5000,
999 int verbose_=1)
1000 : ISTLBackend_OVLP_UMFPack_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
1001 {}
1002 };
1003
1004
1008 template<class GFS>
1010 : public LinearResultStorage
1011 {
1012 public:
1017 explicit ISTLBackend_OVLP_ExplicitDiagonal (const GFS& gfs_)
1018 : gfs(gfs_)
1019 {}
1020
1022 : gfs(other_.gfs)
1023 {}
1024
1029 template<class V>
1030 typename V::ElementType norm(const V& v) const
1031 {
1032 static_assert
1034 "ISTLBackend_OVLP_ExplicitDiagonal::norm() should not be "
1035 "neccessary, so we skipped the implementation. If you have a "
1036 "scenario where you need it, please implement it or report back to "
1037 "us.");
1038 }
1039
1047 template<class M, class V, class W>
1048 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
1049 {
1050 using Backend::Native;
1051 using Backend::native;
1053 Native<M>,
1054 Native<V>,
1055 Native<W>
1056 > jac(native(A),1,1.0);
1057 jac.pre(native(z),native(r));
1058 jac.apply(native(z),native(r));
1059 jac.post(native(z));
1060 if (gfs.gridView().comm().size()>1)
1061 {
1062 CopyDataHandle<GFS,V> copydh(gfs,z);
1063 gfs.gridView().communicate(copydh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
1064 }
1065 res.converged = true;
1066 res.iterations = 1;
1067 res.elapsed = 0.0;
1068 res.reduction = static_cast<double>(reduction);
1069 res.conv_rate = static_cast<double>(reduction); // pow(reduction,1.0/1)
1070 }
1071
1072 private:
1073 const GFS& gfs;
1074 };
1076
1077 template<class GO, int s, template<class,class,class,int> class Preconditioner,
1078 template<class> class Solver>
1079 class ISTLBackend_AMG : public LinearResultStorage
1080 {
1081 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1082 typedef ISTL::ParallelHelper<GFS> PHELPER;
1083 typedef typename GO::Traits::Jacobian M;
1084 typedef Backend::Native<M> MatrixType;
1085 typedef typename GO::Traits::Domain V;
1086 typedef Backend::Native<V> VectorType;
1087 typedef typename ISTL::CommSelector<s,Dune::MPIHelper::isFake>::type Comm;
1088#if HAVE_MPI
1092#else
1095#endif
1096 typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
1098
1099 typedef typename V::ElementType RF;
1100
1101 public:
1102
1107
1108 public:
1109 ISTLBackend_AMG(const GFS& gfs_, unsigned maxiter_=5000,
1110 int verbose_=1, bool reuse_=false,
1111 bool usesuperlu_=true)
1112 : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), params(15,2000),
1113 verbose(verbose_), reuse(reuse_), firstapply(true),
1114 usesuperlu(usesuperlu_)
1115 {
1116 params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
1117 params.setDebugLevel(verbose_);
1118#if !HAVE_SUPERLU
1119 if (gfs.gridView().comm().rank() == 0 && usesuperlu == true)
1120 {
1121 std::cout << "WARNING: You are using AMG without SuperLU!"
1122 << " Please consider installing SuperLU,"
1123 << " or set the usesuperlu flag to false"
1124 << " to suppress this warning." << std::endl;
1125 }
1126#endif
1127 }
1128
1133 void setParameters(const Parameters& params_)
1134 {
1135 params = params_;
1136 }
1137
1145 const Parameters& parameters() const
1146 {
1147 return params;
1148 }
1149
1151 void setReuse(bool reuse_)
1152 {
1153 reuse = reuse_;
1154 }
1155
1157 bool getReuse() const
1158 {
1159 return reuse;
1160 }
1161
1166 typename V::ElementType norm (const V& v) const
1167 {
1168 typedef OverlappingScalarProduct<GFS,V> PSP;
1169 PSP psp(gfs,phelper);
1170 return psp.norm(v);
1171 }
1172
1180 void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
1181 {
1182 Timer watch;
1183 Comm oocc(gfs.gridView().comm());
1184 MatrixType& mat=Backend::native(A);
1186 Dune::Amg::FirstDiagonal> > Criterion;
1187#if HAVE_MPI
1188 phelper.createIndexSetAndProjectForAMG(A, oocc);
1189 Operator oop(mat, oocc);
1191#else
1192 Operator oop(mat);
1194#endif
1195 SmootherArgs smootherArgs;
1196 smootherArgs.iterations = 1;
1197 smootherArgs.relaxationFactor = 1;
1198 Criterion criterion(params);
1199 stats.tprepare=watch.elapsed();
1200 watch.reset();
1201
1202 int verb=0;
1203 if (gfs.gridView().comm().rank()==0) verb=verbose;
1204 //only construct a new AMG if the matrix changes
1205 if (reuse==false || firstapply==true){
1206 amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
1207 firstapply = false;
1208 stats.tsetup = watch.elapsed();
1209 stats.levels = amg->maxlevels();
1210 stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
1211 }
1212 watch.reset();
1213 Solver<VectorType> solver(oop,sp,*amg,RF(reduction),maxiter,verb);
1215
1216 solver.apply(Backend::native(z),Backend::native(r),stat);
1217 stats.tsolve= watch.elapsed();
1218 res.converged = stat.converged;
1219 res.iterations = stat.iterations;
1220 res.elapsed = stat.elapsed;
1221 res.reduction = stat.reduction;
1222 res.conv_rate = stat.conv_rate;
1223 }
1224
1229 const ISTLAMGStatistics& statistics() const
1230 {
1231 return stats;
1232 }
1233
1234 private:
1235 const GFS& gfs;
1236 PHELPER phelper;
1237 unsigned maxiter;
1238 Parameters params;
1239 int verbose;
1240 bool reuse;
1241 bool firstapply;
1242 bool usesuperlu;
1243 std::shared_ptr<AMG> amg;
1244 ISTLAMGStatistics stats;
1245 };
1246
1249
1256 template<class GO, int s=96>
1258 : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1259 {
1260 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1261 public:
1271 ISTLBackend_CG_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1272 int verbose_=1, bool reuse_=false,
1273 bool usesuperlu_=true)
1274 : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1275 (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1276 {}
1277 };
1278
1285 template<class GO, int s=96>
1287 : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1288 {
1289 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1290 public:
1300 ISTLBackend_BCGS_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1301 int verbose_=1, bool reuse_=false,
1302 bool usesuperlu_=true)
1303 : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1304 (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1305 {}
1306 };
1307
1314 template<class GO, int s=96>
1316 : public ISTLBackend_AMG<GO, s, Dune::SeqILU, Dune::BiCGSTABSolver>
1317 {
1318 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1319 public:
1329 ISTLBackend_BCGS_AMG_ILU0(const GFS& gfs_, unsigned maxiter_=5000,
1330 int verbose_=1, bool reuse_=false,
1331 bool usesuperlu_=true)
1332 : ISTLBackend_AMG<GO, s, Dune::SeqILU, Dune::BiCGSTABSolver>
1333 (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1334 {}
1335 };
1336
1338
1340
1341 } // namespace PDELab
1342} // namespace Dune
1343
1344#endif // DUNE_PDELAB_BACKEND_ISTL_OVLPISTLSOLVERBACKEND_HH
The AMG preconditioner.
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:66
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:455
All parameters for AMG.
Definition: parameters.hh:416
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:519
A linear operator exporting itself in matrix form.
Definition: operators.hh:110
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:419
Block parallel preconditioner.
Definition: schwarz.hh:278
conjugate gradient method
Definition: solvers.hh:193
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:135
An overlapping Schwarz operator.
Definition: schwarz.hh:75
Scalar product for overlapping Schwarz methods.
Definition: scalarproducts.hh:201
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by ILU0.
Definition: ovlpistlsolverbackend.hh:1317
ISTLBackend_BCGS_AMG_ILU0(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1329
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR.
Definition: ovlpistlsolverbackend.hh:1288
ISTLBackend_BCGS_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1300
Overlapping parallel conjugate gradient solver preconditioned with AMG smoothed by SSOR.
Definition: ovlpistlsolverbackend.hh:1259
ISTLBackend_CG_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1271
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:685
ISTLBackend_OVLP_BCGS_ILU0(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:694
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:706
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:716
Overlapping parallel BiCGStab solver with SSOR preconditioner.
Definition: ovlpistlsolverbackend.hh:662
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:672
Overlapping parallel BiCGStab solver with SuperLU preconditioner.
Definition: ovlpistlsolverbackend.hh:938
ISTLBackend_OVLP_BCGS_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:948
Overlapping parallel CGS solver with SSOR preconditioner.
Definition: ovlpistlsolverbackend.hh:728
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:738
Overlapping parallel CG solver with SuperLU preconditioner.
Definition: ovlpistlsolverbackend.hh:962
ISTLBackend_OVLP_CG_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:972
Overlapping parallel CG solver with UMFPack preconditioner.
Definition: ovlpistlsolverbackend.hh:987
ISTLBackend_OVLP_CG_UMFPack(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:997
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: ovlpistlsolverbackend.hh:1011
ISTLBackend_OVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:1017
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: ovlpistlsolverbackend.hh:1030
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:1048
Overlapping parallel restarted GMRes solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:752
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:761
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:774
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
virtual void apply(Dune::PDELab::Backend::Vector< GFS, P::domain_type::field_type > &v, const Dune::PDELab::Backend::Vector< GFS, P::range_type::field_type > &d)=0
Apply one step of the preconditioner to the system A(v)=d.
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:827
void apply(X &x, Y &b, InverseOperatorResult &res) override
Apply inverse operator.
Definition: solvers.hh:910
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:52
Sequential ILU preconditioner.
Definition: preconditioners.hh:697
The sequential jacobian preconditioner.
Definition: preconditioners.hh:413
Sequential SSOR preconditioner.
Definition: preconditioners.hh:142
Default implementation for the scalar case.
Definition: scalarproducts.hh:168
Definition: recipe-operator-splitting.cc:108
Definition of the DUNE_NO_DEPRECATED_* macros.
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< IsNumber< A >::value &&!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:42
@ ForwardCommunication
communicate as given in InterfaceType
Definition: gridenums.hh:171
@ InteriorBorder_All_Interface
send interior and border, receive all entities
Definition: gridenums.hh:88
@ All_All_Interface
send all and receive all entities
Definition: gridenums.hh:91
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:13
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:124
The default class for the smoother arguments.
Definition: smoother.hh:38
Statistics about the application of an inverse operator.
Definition: solver.hh:50
double elapsed
Elapsed time in seconds.
Definition: solver.hh:84
int iterations
Number of iterations.
Definition: solver.hh:69
double reduction
Reduction achieved: .
Definition: solver.hh:72
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:78
bool converged
True if convergence criterion has been met.
Definition: solver.hh:75
Category
Definition: solvercategory.hh:23
@ overlapping
Category for overlapping solvers.
Definition: solvercategory.hh:29
Classes for using SuperLU with ISTL matrices.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)