DUNE PDELab (git)

novlpistlsolverbackend.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_NOVLPISTLSOLVERBACKEND_HH
4#define DUNE_PDELAB_BACKEND_ISTL_NOVLPISTLSOLVERBACKEND_HH
5
6#include <cstddef>
7
10
11#include <dune/grid/common/gridenums.hh>
12
13#include <dune/istl/io.hh>
17#include <dune/istl/paamg/pinfo.hh>
20#include <dune/istl/solvercategory.hh>
21#include <dune/istl/solvers.hh>
22#include <dune/istl/superlu.hh>
23
24#include <dune/pdelab/constraints/common/constraints.hh>
25#include <dune/pdelab/gridfunctionspace/genericdatahandle.hh>
26#include <dune/pdelab/backend/istl/vector.hh>
27#include <dune/pdelab/backend/istl/bcrsmatrix.hh>
28#include <dune/pdelab/backend/istl/blockmatrixdiagonal.hh>
29#include <dune/pdelab/backend/istl/parallelhelper.hh>
30#include <dune/pdelab/backend/istl/seqistlsolverbackend.hh>
31
32namespace Dune {
33 namespace PDELab {
34
38
39 //========================================================
40 // Generic support for nonoverlapping grids
41 //========================================================
42
44
52 template<typename GFS, typename M, typename X, typename Y>
54 : public Dune::AssembledLinearOperator<M,X,Y>
55 {
56 public:
58 using matrix_type = Backend::Native<M>;
60 using domain_type = Backend::Native<X>;
62 using range_type = Backend::Native<Y>;
64 typedef typename X::field_type field_type;
65
67
77 NonoverlappingOperator (const GFS& gfs_, const M& A)
78 : gfs(gfs_), _A_(A)
79 { }
80
82
87 virtual void apply (const X& x, Y& y) const override
88 {
89 using Backend::native;
90 // apply local operator; now we have sum y_p = sequential y
91 native(_A_).mv(native(x),native(y));
92
93 // accumulate y on border
94 Dune::PDELab::AddDataHandle<GFS,Y> adddh(gfs,y);
95 if (gfs.gridView().comm().size()>1)
97 }
98
100
105 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
106 {
107 using Backend::native;
108 // apply local operator; now we have sum y_p = sequential y
109 native(_A_).usmv(alpha,native(x),native(y));
110
111 // accumulate y on border
112 Dune::PDELab::AddDataHandle<GFS,Y> adddh(gfs,y);
113 if (gfs.gridView().comm().size()>1)
115 }
116
118 {
120 }
121
123 virtual const M& getmat () const override
124 {
125 return _A_;
126 }
127
128 private:
129 const GFS& gfs;
130 const M& _A_;
131 };
132
133 // parallel scalar product assuming no overlap
134 template<class GFS, class X>
135 class NonoverlappingScalarProduct : public Dune::ScalarProduct<X>
136 {
137 public:
139 typedef X domain_type;
140 typedef typename X::ElementType field_type;
141
142 SolverCategory::Category category() const override
143 {
145 }
146
149 NonoverlappingScalarProduct (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_)
150 : gfs(gfs_), helper(helper_)
151 {}
152
157 virtual field_type dot (const X& x, const X& y) const override
158 {
159 // do local scalar product on unique partition
160 field_type sum = helper.disjointDot(x,y);
161
162 // do global communication
163 return gfs.gridView().comm().sum(sum);
164 }
165
169 virtual double norm (const X& x) const override
170 {
171 using std::sqrt;
172 return sqrt(static_cast<double>(this->dot(x,x)));
173 }
174
177 void make_consistent (X& x) const
178 {
179 Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,x);
180 if (gfs.gridView().comm().size()>1)
182 }
183
184 private:
185 const GFS& gfs;
186 const ISTL::ParallelHelper<GFS>& helper;
187 };
188
189 // parallel Richardson preconditioner
190 template<class GFS, class X, class Y>
191 class NonoverlappingRichardson : public Dune::Preconditioner<X,Y>
192 {
193 public:
195 typedef X domain_type;
197 typedef Y range_type;
199 typedef typename X::ElementType field_type;
200
201 // define the category
202 SolverCategory::Category category() const override
203 {
205 }
206
208 NonoverlappingRichardson (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_)
209 : gfs(gfs_), helper(helper_)
210 {
211 }
212
216 virtual void pre (X& x, Y& b) const override {}
217
221 virtual void apply (X& v, const Y& d) const override
222 {
223 v = d;
224 }
225
229 virtual void post (X& x) override {}
230
231 private:
232 const GFS& gfs;
233 const ISTL::ParallelHelper<GFS>& helper;
234 };
235
237
250 template<typename A, typename X, typename Y>
252 : public Dune::Preconditioner<X,Y>
253 {
254
255 typedef typename ISTL::BlockMatrixDiagonal<A>::MatrixElementVector Diagonal;
256
257 Diagonal _inverse_diagonal;
258
259 public:
261
265 typedef X domain_type;
267
271 typedef Y range_type;
273 typedef typename X::ElementType field_type;
274
276 {
278 }
279
281
292 template<typename GFS>
293 NonoverlappingJacobi(const GFS& gfs, const A &m)
294 : _inverse_diagonal(m)
295 {
296 // make the diagonal consistent...
297 typename ISTL::BlockMatrixDiagonal<A>::template AddMatrixElementVectorDataHandle<GFS> addDH(gfs, _inverse_diagonal);
298 gfs.gridView().communicate(addDH,
301
302 // ... and then invert it
303 _inverse_diagonal.invert();
304
305 }
306
308 virtual void pre (X& x, Y& b) override {}
309
311 /*
312 * For this preconditioner, this method works with both consistent and
313 * inconsistent vectors: if d is consistent, v will be consistent, if d
314 * is inconsistent, v will be inconsistent.
315 */
316 virtual void apply (X& v, const Y& d) override
317 {
318 _inverse_diagonal.mv(d,v);
319 }
320
322 virtual void post (X& x) override {}
323 };
324
327
329 template<class GFS>
331 {
332 typedef ISTL::ParallelHelper<GFS> PHELPER;
333
334 public:
341 explicit ISTLBackend_NOVLP_CG_NOPREC (const GFS& gfs_,
342 unsigned maxiter_=5000,
343 int verbose_=1)
344 : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
345 {}
346
351 template<class V>
352 typename V::ElementType norm (const V& v) const
353 {
354 V x(v); // make a copy because it has to be made consistent
355 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
356 PSP psp(gfs,phelper);
357 psp.make_consistent(x);
358 return psp.norm(x);
359 }
360
368 template<class M, class V, class W>
369 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
370 {
372 POP pop(gfs,A);
373 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
374 PSP psp(gfs,phelper);
375 typedef Dune::PDELab::NonoverlappingRichardson<GFS,V,W> PRICH;
376 PRICH prich(gfs,phelper);
377 int verb=0;
378 if (gfs.gridView().comm().rank()==0) verb=verbose;
379 Dune::CGSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
381 solver.apply(z,r,stat);
382 res.converged = stat.converged;
383 res.iterations = stat.iterations;
384 res.elapsed = stat.elapsed;
385 res.reduction = stat.reduction;
386 res.conv_rate = stat.conv_rate;
387 }
388
390 const Dune::PDELab::LinearSolverResult<double>& result() const
391 {
392 return res;
393 }
394
395 private:
396 const GFS& gfs;
397 PHELPER phelper;
398 Dune::PDELab::LinearSolverResult<double> res;
399 unsigned maxiter;
400 int verbose;
401 };
402
404 template<class GFS>
406 {
407 typedef ISTL::ParallelHelper<GFS> PHELPER;
408
409 const GFS& gfs;
410 PHELPER phelper;
411 LinearSolverResult<double> res;
412 unsigned maxiter;
413 int verbose;
414
415 public:
417
422 explicit ISTLBackend_NOVLP_CG_Jacobi(const GFS& gfs_,
423 unsigned maxiter_ = 5000,
424 int verbose_ = 1) :
425 gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
426 {}
427
429
434 template<class V>
435 typename V::ElementType norm (const V& v) const
436 {
437 V x(v); // make a copy because it has to be made consistent
438 typedef NonoverlappingScalarProduct<GFS,V> PSP;
439 PSP psp(gfs,phelper);
440 psp.make_consistent(x);
441 return psp.norm(x);
442 }
443
445
457 template<class M, class V, class W>
458 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
459 {
461 POP pop(gfs,A);
462 typedef NonoverlappingScalarProduct<GFS,V> PSP;
463 PSP psp(gfs,phelper);
464
465 typedef NonoverlappingJacobi<M,V,W> PPre;
466 PPre ppre(gfs,Backend::native(A));
467
468 int verb=0;
469 if (gfs.gridView().comm().rank()==0) verb=verbose;
470 CGSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
472 solver.apply(z,r,stat);
473 res.converged = stat.converged;
474 res.iterations = stat.iterations;
475 res.elapsed = stat.elapsed;
476 res.reduction = stat.reduction;
477 res.conv_rate = stat.conv_rate;
478 }
479
481 const LinearSolverResult<double>& result() const
482 { return res; }
483 };
484
486 template<class GFS>
488 {
489 typedef ISTL::ParallelHelper<GFS> PHELPER;
490
491 public:
498 explicit ISTLBackend_NOVLP_BCGS_NOPREC (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
499 : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
500 {}
501
506 template<class V>
507 typename V::ElementType norm (const V& v) const
508 {
509 V x(v); // make a copy because it has to be made consistent
510 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
511 PSP psp(gfs,phelper);
512 psp.make_consistent(x);
513 return psp.norm(x);
514 }
515
523 template<class M, class V, class W>
524 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
525 {
527 POP pop(gfs,A);
528 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
529 PSP psp(gfs,phelper);
530 typedef Dune::PDELab::NonoverlappingRichardson<GFS,V,W> PRICH;
531 PRICH prich(gfs,phelper);
532 int verb=0;
533 if (gfs.gridView().comm().rank()==0) verb=verbose;
534 Dune::BiCGSTABSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
536 solver.apply(z,r,stat);
537 res.converged = stat.converged;
538 res.iterations = stat.iterations;
539 res.elapsed = stat.elapsed;
540 res.reduction = stat.reduction;
541 res.conv_rate = stat.conv_rate;
542 }
543
545 const Dune::PDELab::LinearSolverResult<double>& result() const
546 {
547 return res;
548 }
549
550 private:
551 const GFS& gfs;
552 PHELPER phelper;
553 Dune::PDELab::LinearSolverResult<double> res;
554 unsigned maxiter;
555 int verbose;
556 };
557
558
560 template<class GFS>
562 {
563 typedef ISTL::ParallelHelper<GFS> PHELPER;
564
565 public:
572 explicit ISTLBackend_NOVLP_BCGS_Jacobi (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
573 : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
574 {}
575
580 template<class V>
581 typename V::ElementType norm (const V& v) const
582 {
583 V x(v); // make a copy because it has to be made consistent
584 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
585 PSP psp(gfs,phelper);
586 psp.make_consistent(x);
587 return psp.norm(x);
588 }
589
597 template<class M, class V, class W>
598 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
599 {
601 POP pop(gfs,A);
602 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
603 PSP psp(gfs,phelper);
604
605 typedef NonoverlappingJacobi<M,V,W> PPre;
606 PPre ppre(gfs,A);
607
608 int verb=0;
609 if (gfs.gridView().comm().rank()==0) verb=verbose;
610 Dune::BiCGSTABSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
612 solver.apply(z,r,stat);
613 res.converged = stat.converged;
614 res.iterations = stat.iterations;
615 res.elapsed = stat.elapsed;
616 res.reduction = stat.reduction;
617 res.conv_rate = stat.conv_rate;
618 }
619
621 const Dune::PDELab::LinearSolverResult<double>& result() const
622 {
623 return res;
624 }
625
626 private:
627 const GFS& gfs;
628 PHELPER phelper;
629 Dune::PDELab::LinearSolverResult<double> res;
630 unsigned maxiter;
631 int verbose;
632 };
633
635 template<typename GFS>
637 {
638 typedef ISTL::ParallelHelper<GFS> PHELPER;
639
640 const GFS& gfs;
641 PHELPER phelper;
642 Dune::PDELab::LinearSolverResult<double> res;
643
644 public:
650 explicit ISTLBackend_NOVLP_ExplicitDiagonal(const GFS& gfs_)
651 : gfs(gfs_), phelper(gfs)
652 {}
653
658 template<class V>
659 typename V::ElementType norm (const V& v) const
660 {
661 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
662 V x(v); // make a copy because it has to be made consistent
663 PSP psp(gfs,phelper);
664 psp.make_consistent(x);
665 return psp.norm(x);
666 }
667
675 template<class M, class V, class W>
676 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
677 {
678 Dune::SeqJac<M,V,W> jac(A,1,1.0);
679 jac.pre(z,r);
680 jac.apply(z,r);
681 jac.post(z);
682 if (gfs.gridView().comm().size()>1)
683 {
684 Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,z);
686 }
687 res.converged = true;
688 res.iterations = 1;
689 res.elapsed = 0.0;
690 res.reduction = reduction;
691 res.conv_rate = reduction; // pow(reduction,1.0/1)
692 }
693
695 const Dune::PDELab::LinearSolverResult<double>& result() const
696 {
697 return res;
698 }
699 };
701
702
710 template<class GO,
711 template<class,class,class,int> class Preconditioner,
712 template<class> class Solver>
714 {
715 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
716 typedef ISTL::ParallelHelper<GFS> PHELPER;
717
718 public:
726 explicit ISTLBackend_NOVLP_BASE_PREC (const GO& grid_operator, unsigned maxiter_ = 5000, unsigned steps_ = 5, int verbose_ = 1)
727 : _grid_operator(grid_operator)
728 , gfs(grid_operator.trialGridFunctionSpace())
729 , phelper(gfs,verbose_)
730 , maxiter(maxiter_)
731 , steps(steps_)
732 , verbose(verbose_)
733 {}
734
739 template<class Vector>
740 typename Vector::ElementType norm (const Vector& v) const
741 {
742 Vector x(v); // make a copy because it has to be made consistent
743 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,Vector> PSP;
744 PSP psp(gfs,phelper);
745 psp.make_consistent(x);
746 return psp.norm(x);
747 }
748
756 template<class M, class V, class W>
757 void apply(M& A, V& z, W& r, typename V::ElementType reduction)
758 {
759 using MatrixType = Backend::Native<M>;
760 MatrixType& mat = Backend::native(A);
761 using VectorType = Backend::Native<W>;
762#if HAVE_MPI
763 typedef typename ISTL::CommSelector<96,Dune::MPIHelper::isFake>::type Comm;
764 _grid_operator.make_consistent(A);
765 ISTL::assertParallelUG(gfs.gridView().comm());
766 Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
767 phelper.createIndexSetAndProjectForAMG(mat, oocc);
769 Smoother smoother(mat, steps, 1.0);
771 PSP psp(oocc);
773 Operator oop(mat,oocc);
775 ParSmoother parsmoother(smoother, oocc);
776#else
778 ParSmoother parsmoother(mat, steps, 1.0);
780 PSP psp;
782 Operator oop(mat);
783#endif
784 int verb=0;
785 if (gfs.gridView().comm().rank()==0) verb=verbose;
786 Solver<VectorType> solver(oop,psp,parsmoother,reduction,maxiter,verb);
788 //make r consistent
789 if (gfs.gridView().comm().size()>1){
790 Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,r);
791 gfs.gridView().communicate(adddh,
794 }
795
796 solver.apply(z,r,stat);
797 res.converged = stat.converged;
798 res.iterations = stat.iterations;
799 res.elapsed = stat.elapsed;
800 res.reduction = stat.reduction;
801 res.conv_rate = stat.conv_rate;
802 }
803
805 const Dune::PDELab::LinearSolverResult<double>& result() const
806 {
807 return res;
808 }
809
810 private:
811 const GO& _grid_operator;
812 const GFS& gfs;
813 PHELPER phelper;
814 Dune::PDELab::LinearSolverResult<double> res;
815 unsigned maxiter;
816 unsigned steps;
817 int verbose;
818 };
819
822
835 template<class GO>
837 : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>
838 {
839
840 public:
848 explicit ISTLBackend_NOVLP_BCGS_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
849 int steps_=5, int verbose_=1)
850 : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_, steps_, verbose_)
851 {}
852 };
853
860 template<class GO>
862 : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>
863 {
864
865 public:
873 explicit ISTLBackend_NOVLP_CG_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
874 int steps_=5, int verbose_=1)
875 : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_, steps_, verbose_)
876 {}
877 };
880
881 template<class GO,int s, template<class,class,class,int> class Preconditioner,
882 template<class> class Solver>
883 class ISTLBackend_AMG_NOVLP : public LinearResultStorage
884 {
885 typedef typename GO::Traits::TrialGridFunctionSpace GFS;
886 typedef typename ISTL::ParallelHelper<GFS> PHELPER;
887 typedef typename GO::Traits::Jacobian M;
888 typedef Backend::Native<M> MatrixType;
889 typedef typename GO::Traits::Domain V;
890 typedef Backend::Native<V> VectorType;
891 typedef typename ISTL::CommSelector<s,Dune::MPIHelper::isFake>::type Comm;
892#if HAVE_MPI
896#else
899#endif
900 typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
903
904 public:
905 ISTLBackend_AMG_NOVLP(const GO& grid_operator, unsigned maxiter_=5000,
906 int verbose_=1, bool reuse_=false,
907 bool usesuperlu_=true)
908 : _grid_operator(grid_operator)
909 , gfs(grid_operator.trialGridFunctionSpace())
910 , phelper(gfs,verbose_)
911 , maxiter(maxiter_)
912 , params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu)
913 , verbose(verbose_)
914 , reuse(reuse_)
915 , firstapply(true)
916 , usesuperlu(usesuperlu_)
917 {
918 params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
919 params.setDebugLevel(verbose_);
920#if !HAVE_SUPERLU
921 if (phelper.rank() == 0 && usesuperlu == true)
922 {
923 std::cout << "WARNING: You are using AMG without SuperLU!"
924 << " Please consider installing SuperLU,"
925 << " or set the usesuperlu flag to false"
926 << " to suppress this warning." << std::endl;
927 }
928#endif
929 }
930
935 void setParameters(const Parameters& params_)
936 {
937 params = params_;
938 }
939
947 const Parameters& parameters() const
948 {
949 return params;
950 }
951
953 void setReuse(bool reuse_)
954 {
955 reuse = reuse_;
956 }
957
959 bool getReuse() const
960 {
961 return reuse;
962 }
963
964
969 typename V::ElementType norm (const V& v) const
970 {
971 V x(v); // make a copy because it has to be made consistent
972 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
973 PSP psp(gfs,phelper);
974 psp.make_consistent(x);
975 return psp.norm(x);
976 }
977
978 void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
979 {
980 Timer watch;
981 MatrixType& mat = Backend::native(A);
983 Dune::Amg::FirstDiagonal> > Criterion;
984#if HAVE_MPI
985 Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
986 _grid_operator.make_consistent(A);
987 phelper.createIndexSetAndProjectForAMG(A, oocc);
989 Operator oop(mat, oocc);
990#else
991 Comm oocc(gfs.gridView().comm());
992 Operator oop(mat);
994#endif
995 SmootherArgs smootherArgs;
996 smootherArgs.iterations = 1;
997 smootherArgs.relaxationFactor = 1;
998 //use noAccu or atOnceAccu
999 Criterion criterion(params);
1000 stats.tprepare=watch.elapsed();
1001 watch.reset();
1002
1003 int verb=0;
1004 if (gfs.gridView().comm().rank()==0) verb=verbose;
1005 //only construct a new AMG if the matrix changes
1006 if (reuse==false || firstapply==true){
1007 amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
1008 firstapply = false;
1009 stats.tsetup = watch.elapsed();
1010 stats.levels = amg->maxlevels();
1011 stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
1012 }
1013
1015 // make r consistent
1016 if (gfs.gridView().comm().size()>1) {
1017 Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,r);
1018 gfs.gridView().communicate(adddh,
1021 }
1022 watch.reset();
1023 Solver<VectorType> solver(oop,sp,*amg,reduction,maxiter,verb);
1024 solver.apply(Backend::native(z),Backend::native(r),stat);
1025 stats.tsolve= watch.elapsed();
1026 res.converged = stat.converged;
1027 res.iterations = stat.iterations;
1028 res.elapsed = stat.elapsed;
1029 res.reduction = stat.reduction;
1030 res.conv_rate = stat.conv_rate;
1031 }
1032
1037 const ISTLAMGStatistics& statistics() const
1038 {
1039 return stats;
1040 }
1041
1042 private:
1043 const GO& _grid_operator;
1044 const GFS& gfs;
1045 PHELPER phelper;
1046 unsigned maxiter;
1047 Parameters params;
1048 int verbose;
1049 bool reuse;
1050 bool firstapply;
1051 bool usesuperlu;
1052 std::shared_ptr<AMG> amg;
1053 ISTLAMGStatistics stats;
1054 };
1055
1058
1073 template<class GO, int s=96>
1075 : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1076 {
1077
1078 public:
1079 ISTLBackend_NOVLP_CG_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1080 int verbose_=1, bool reuse_=false,
1081 bool usesuperlu_=true)
1082 : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1083 {}
1084 };
1085
1100 template<class GO, int s=96>
1102 : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1103 {
1104
1105 public:
1106 ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1107 int verbose_=1, bool reuse_=false,
1108 bool usesuperlu_=true)
1109 : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1110 {}
1111 };
1112
1127 template<class GO, int s=96>
1129 : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>
1130 {
1131
1132 public:
1133 ISTLBackend_NOVLP_LS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1134 int verbose_=1, bool reuse_=false,
1135 bool usesuperlu_=true)
1136 : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1137 {}
1138 };
1141
1142 } // namespace PDELab
1143} // namespace Dune
1144
1145#endif // DUNE_PDELAB_BACKEND_ISTL_NOVLPISTLSOLVERBACKEND_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
M matrix_type
export types, usually they come from the derived class
Definition: operators.hh:113
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:419
void apply(X &x, X &b, InverseOperatorResult &res) override
Apply inverse operator.
Definition: solvers.hh:439
conjugate gradient method
Definition: solvers.hh:193
void apply(X &x, X &b, InverseOperatorResult &res) override
Apply inverse operator.
Definition: solvers.hh:279
X::field_type field_type
The field type of the operator.
Definition: operators.hh:75
Y range_type
The type of the range of the operator.
Definition: operators.hh:73
X domain_type
The type of the domain of the operator.
Definition: operators.hh:71
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:135
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:276
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:61
Nonoverlapping Scalar Product with communication object.
Definition: scalarproducts.hh:179
Utility base class for preconditioned novlp backends.
Definition: novlpistlsolverbackend.hh:714
ISTLBackend_NOVLP_BASE_PREC(const GO &grid_operator, unsigned maxiter_=5000, unsigned steps_=5, int verbose_=1)
Constructor.
Definition: novlpistlsolverbackend.hh:726
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:805
Vector::ElementType norm(const Vector &v) const
Compute global norm of a vector.
Definition: novlpistlsolverbackend.hh:740
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
Solve the given linear system.
Definition: novlpistlsolverbackend.hh:757
Nonoverlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1103
Nonoverlapping parallel BiCGStab solver with Jacobi preconditioner.
Definition: novlpistlsolverbackend.hh:562
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:621
ISTLBackend_NOVLP_BCGS_Jacobi(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:572
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:581
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:598
Nonoverlapping parallel BiCGStab solver without preconditioner.
Definition: novlpistlsolverbackend.hh:488
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:545
ISTLBackend_NOVLP_BCGS_NOPREC(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:498
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:524
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:507
Nonoverlapping parallel BiCGSTAB solver preconditioned by block SSOR.
Definition: novlpistlsolverbackend.hh:838
ISTLBackend_NOVLP_BCGS_SSORk(const GO &grid_operator, unsigned maxiter_=5000, int steps_=5, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:848
Nonoverlapping parallel CG solver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1076
Nonoverlapping parallel CG solver with Jacobi preconditioner.
Definition: novlpistlsolverbackend.hh:406
const LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:481
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:435
ISTLBackend_NOVLP_CG_Jacobi(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:422
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:458
Nonoverlapping parallel CG solver without preconditioner.
Definition: novlpistlsolverbackend.hh:331
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:369
ISTLBackend_NOVLP_CG_NOPREC(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:341
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:390
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:352
Nonoverlapping parallel CG solver preconditioned by block SSOR.
Definition: novlpistlsolverbackend.hh:863
ISTLBackend_NOVLP_CG_SSORk(const GO &grid_operator, unsigned maxiter_=5000, int steps_=5, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:873
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: novlpistlsolverbackend.hh:637
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:659
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:695
ISTLBackend_NOVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: novlpistlsolverbackend.hh:650
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename W::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:676
Nonoverlapping parallel LoopSolver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1130
parallel non-overlapping Jacobi preconditioner
Definition: novlpistlsolverbackend.hh:253
virtual void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition: novlpistlsolverbackend.hh:308
virtual void post(X &x) override
Clean up.
Definition: novlpistlsolverbackend.hh:322
virtual void apply(X &v, const Y &d) override
Apply the precondioner.
Definition: novlpistlsolverbackend.hh:316
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
Definition: novlpistlsolverbackend.hh:275
X::ElementType field_type
The field type of the preconditioner.
Definition: novlpistlsolverbackend.hh:273
X domain_type
The domain type of the operator.
Definition: novlpistlsolverbackend.hh:265
NonoverlappingJacobi(const GFS &gfs, const A &m)
Constructor.
Definition: novlpistlsolverbackend.hh:293
Y range_type
The range type of the operator.
Definition: novlpistlsolverbackend.hh:271
Operator for the non-overlapping parallel case.
Definition: novlpistlsolverbackend.hh:55
NonoverlappingOperator(const GFS &gfs_, const M &A)
Construct a non-overlapping operator.
Definition: novlpistlsolverbackend.hh:77
virtual const M & getmat() const override
extract the matrix
Definition: novlpistlsolverbackend.hh:123
virtual void apply(const X &x, Y &y) const override
apply operator
Definition: novlpistlsolverbackend.hh:87
X::field_type field_type
export type of the entries for x
Definition: novlpistlsolverbackend.hh:64
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const override
apply operator to x, scale and add:
Definition: novlpistlsolverbackend.hh:105
SolverCategory::Category category() const override
Category of the linear operator (see SolverCategory::Category)
Definition: novlpistlsolverbackend.hh:117
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
virtual void apply(X &v, const Y &d)=0
Apply one step of the preconditioner to the system A(v)=d.
virtual void pre(X &x, Y &b)=0
Prepare the preconditioner.
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:52
The sequential jacobian preconditioner.
Definition: preconditioners.hh:413
void apply(X &v, const Y &d) override
Apply the preconditioner.
Definition: preconditioners.hh:488
void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition: preconditioners.hh:480
void post(X &x) override
Clean up.
Definition: preconditioners.hh:500
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_InteriorBorder_Interface
send/receive interior and border entities
Definition: gridenums.hh:87
@ atOnceAccu
Accumulate data to one process at once.
Definition: parameters.hh:243
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.
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
@ nonoverlapping
Category for non-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 (Nov 23, 23:29, 2024)