dune-pdelab  2.4.1
istl/ovlpistlsolverbackend.hh
Go to the documentation of this file.
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_OVLPISTLSOLVERBACKEND_HH
4 #define DUNE_OVLPISTLSOLVERBACKEND_HH
5 
6 #include <dune/common/deprecated.hh>
7 #include <dune/common/parallel/mpihelper.hh>
8 
9 #include <dune/istl/owneroverlapcopy.hh>
10 #include <dune/istl/solvercategory.hh>
11 #include <dune/istl/operators.hh>
12 #include <dune/istl/solvers.hh>
13 #include <dune/istl/preconditioners.hh>
14 #include <dune/istl/scalarproducts.hh>
15 #include <dune/istl/paamg/amg.hh>
16 #include <dune/istl/paamg/pinfo.hh>
17 #include <dune/istl/io.hh>
18 #include <dune/istl/superlu.hh>
19 
26 
27 namespace 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>
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  //redefine the category, that is the only difference
52  enum {category=Dune::SolverCategory::overlapping};
53 
54  OverlappingOperator (const CC& cc_, const M& A)
55  : cc(cc_), _A_(A)
56  {}
57 
59  virtual void apply (const domain_type& x, range_type& y) const
60  {
61  using Backend::native;
62  native(_A_).mv(native(x),native(y));
64  }
65 
67  virtual void applyscaleadd (field_type alpha, const domain_type& x, range_type& y) const
68  {
69  using Backend::native;
70  native(_A_).usmv(alpha,native(x),native(y));
72  }
73 
75  virtual const M& getmat () const
76  {
77  return _A_;
78  }
79 
80  private:
81  const CC& cc;
82  const M& _A_;
83  };
84 
85  // new scalar product assuming at least overlap 1
86  // uses unique partitioning of nodes for parallelization
87  template<class GFS, class X>
89  : public Dune::ScalarProduct<X>
90  {
91  public:
93  typedef X domain_type;
94  typedef typename X::ElementType field_type;
95 
97  enum {category=Dune::SolverCategory::overlapping};
98 
101  OverlappingScalarProduct (const GFS& gfs_, const istl::ParallelHelper<GFS>& helper_)
102  : gfs(gfs_), helper(helper_)
103  {}
104 
105 
110  virtual field_type dot (const X& x, const X& y)
111  {
112  // do local scalar product on unique partition
113  field_type sum = helper.disjointDot(x,y);
114 
115  // do global communication
116  return gfs.gridView().comm().sum(sum);
117  }
118 
122  virtual double norm (const X& x)
123  {
124  return sqrt(static_cast<double>(this->dot(x,x)));
125  }
126 
127  private:
128  const GFS& gfs;
129  const istl::ParallelHelper<GFS>& helper;
130  };
131 
132  // wrapped sequential preconditioner
133  template<class CC, class GFS, class P>
135  : public Dune::Preconditioner<Dune::PDELab::Backend::Vector<GFS,typename P::domain_type::field_type>,
136  Dune::PDELab::Backend::Vector<GFS,typename P::range_type::field_type>>
137  {
138  public:
143 
144  // define the category
145  enum {
147  category=Dune::SolverCategory::overlapping
148  };
149 
151  OverlappingWrappedPreconditioner (const GFS& gfs_, P& prec_, const CC& cc_,
152  const istl::ParallelHelper<GFS>& helper_)
153  : gfs(gfs_), prec(prec_), cc(cc_), helper(helper_)
154  {}
155 
159  virtual void pre (domain_type& x, range_type& b)
160  {
161  prec.pre(x,b);
162  }
163 
167  virtual void apply (domain_type& v, const range_type& d)
168  {
169  range_type dd(d);
170  set_constrained_dofs(cc,0.0,dd);
171  prec.apply(Backend::native(v),Backend::native(dd));
173  if (gfs.gridView().comm().size()>1)
174  gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
175  }
176 
180  virtual void post (domain_type& x)
181  {
182  prec.post(Backend::native(x));
183  }
184 
185  private:
186  const GFS& gfs;
187  P& prec;
188  const CC& cc;
189  const istl::ParallelHelper<GFS>& helper;
190  };
191 
192 
193 #if HAVE_SUPERLU
194  // exact subdomain solves with SuperLU as preconditioner
195  template<class GFS, class M, class X, class Y>
196  class SuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
197  {
198  typedef Backend::Native<M> ISTLM;
199 
200  public:
202  typedef X domain_type;
204  typedef Y range_type;
206  typedef typename X::ElementType field_type;
207 
208 
209  // define the category
210  enum {
212  category=Dune::SolverCategory::overlapping
213  };
214 
221  SuperLUSubdomainSolver (const GFS& gfs_, const M& A_)
222  : gfs(gfs_), solver(Backend::native(A_),false) // this does the decomposition
223  {}
224 
228  virtual void pre (X& x, Y& b) {}
229 
233  virtual void apply (X& v, const Y& d)
234  {
235  Dune::InverseOperatorResult stat;
236  Y b(d); // need copy, since solver overwrites right hand side
237  solver.apply(Backend::native(v),Backend::native(b),stat);
238  if (gfs.gridView().comm().size()>1)
239  {
240  AddDataHandle<GFS,X> adddh(gfs,v);
241  gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
242  }
243  }
244 
248  virtual void post (X& x) {}
249 
250  private:
251  const GFS& gfs;
252  Dune::SuperLU<ISTLM> solver;
253  };
254 
255  // exact subdomain solves with SuperLU as preconditioner
256  template<class GFS, class M, class X, class Y>
257  class RestrictedSuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
258  {
259  typedef typename M::BaseT ISTLM;
260 
261  public:
263  typedef X domain_type;
265  typedef Y range_type;
267  typedef typename X::ElementType field_type;
268 
269 
270  // define the category
271  enum {
273  category=Dune::SolverCategory::overlapping
274  };
275 
283  RestrictedSuperLUSubdomainSolver (const GFS& gfs_, const M& A_,
284  const istl::ParallelHelper<GFS>& helper_)
285  : gfs(gfs_), solver(Backend::native(A_),false), helper(helper_) // this does the decomposition
286  {}
287 
291  virtual void pre (X& x, Y& b) {}
292 
296  virtual void apply (X& v, const Y& d)
297  {
298  using Backend::native;
299  Dune::InverseOperatorResult stat;
300  Y b(d); // need copy, since solver overwrites right hand side
301  solver.apply(native(v),native(b),stat);
302  if (gfs.gridView().comm().size()>1)
303  {
304  helper.maskForeignDOFs(native(v));
305  AddDataHandle<GFS,X> adddh(gfs,v);
306  gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
307  }
308  }
309 
313  virtual void post (X& x) {}
314 
315  private:
316  const GFS& gfs;
317  Dune::SuperLU<ISTLM> solver;
318  const istl::ParallelHelper<GFS>& helper;
319  };
320 #endif
321 
322  template<typename GFS>
324  {
325  public:
327  : gfs(gfs_), helper(gfs_)
328  {}
329 
334  template<typename X>
335  typename X::ElementType dot (const X& x, const X& y) const
336  {
337  // do local scalar product on unique partition
338  typename X::ElementType sum = helper.disjointDot(x,y);
339 
340  // do global communication
341  return gfs.gridView().comm().sum(sum);
342  }
343 
347  template<typename X>
348  typename Dune::template FieldTraits<typename X::ElementType >::real_type norm (const X& x) const
349  {
350  using namespace std;
351  return sqrt(static_cast<double>(this->dot(x,x)));
352  }
353 
355  {
356  return helper;
357  }
358 
359  // need also non-const version;
360  istl::ParallelHelper<GFS>& parallelHelper() // P.B.: needed for createIndexSetAndProjectForAMG
361  {
362  return helper;
363  }
364 
365  private:
366  const GFS& gfs;
368  };
369 
370 
371  template<typename GFS, typename X>
373  : public ScalarProduct<X>
374  {
375  public:
376  enum {category=Dune::SolverCategory::overlapping};
378  : implementation(implementation_)
379  {}
380 
381  virtual typename X::BaseT::field_type dot(const X& x, const X& y)
382  {
383  return implementation.dot(x,y);
384  }
385 
386  virtual typename X::BaseT::field_type norm (const X& x)
387  {
388  using namespace std;
389  return sqrt(static_cast<double>(this->dot(x,x)));
390  }
391 
392  private:
393  const OVLPScalarProductImplementation<GFS>& implementation;
394  };
395 
396  template<class GFS, class C,
397  template<class,class,class,int> class Preconditioner,
398  template<class> class Solver>
401  {
402  public:
411  ISTLBackend_OVLP_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
412  int steps_=5, int verbose_=1)
413  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), steps(steps_), verbose(verbose_)
414  {}
415 
423  template<class M, class V, class W>
424  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
425  {
426  using Backend::Native;
427  using Backend::native;
428  typedef OverlappingOperator<C,M,V,W> POP;
429  POP pop(c,A);
430  typedef OVLPScalarProduct<GFS,V> PSP;
431  PSP psp(*this);
432  typedef Preconditioner<
433  Native<M>,
434  Native<V>,
435  Native<W>,
436  1
437  > SeqPrec;
438  SeqPrec seqprec(native(A),steps,1.0);
440  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
441  int verb=0;
442  if (gfs.gridView().comm().rank()==0) verb=verbose;
443  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
444  Dune::InverseOperatorResult stat;
445  solver.apply(z,r,stat);
446  res.converged = stat.converged;
447  res.iterations = stat.iterations;
448  res.elapsed = stat.elapsed;
449  res.reduction = stat.reduction;
450  res.conv_rate = stat.conv_rate;
451  }
452  private:
453  const GFS& gfs;
454  const C& c;
455  unsigned maxiter;
456  int steps;
457  int verbose;
458  };
459 
460  // Base class for ILU0 as preconditioner
461  template<class GFS, class C,
462  template<class> class Solver>
465  {
466  public:
474  ISTLBackend_OVLP_ILU0_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000, int verbose_=1)
475  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
476  {}
477 
485  template<class M, class V, class W>
486  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
487  {
488  using Backend::Native;
489  using Backend::native;
490  typedef OverlappingOperator<C,M,V,W> POP;
491  POP pop(c,A);
492  typedef OVLPScalarProduct<GFS,V> PSP;
493  PSP psp(*this);
494  typedef SeqILU0<
495  Native<M>,
496  Native<V>,
497  Native<W>,
498  1
499  > SeqPrec;
500  SeqPrec seqprec(native(A),1.0);
502  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
503  int verb=0;
504  if (gfs.gridView().comm().rank()==0) verb=verbose;
505  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
506  Dune::InverseOperatorResult stat;
507  solver.apply(z,r,stat);
508  res.converged = stat.converged;
509  res.iterations = stat.iterations;
510  res.elapsed = stat.elapsed;
511  res.reduction = stat.reduction;
512  res.conv_rate = stat.conv_rate;
513  }
514  private:
515  const GFS& gfs;
516  const C& c;
517  unsigned maxiter;
518  int steps;
519  int verbose;
520  };
521 
522  // Base class for ILUn as preconditioner
523  template<class GFS, class C,
524  template<class> class Solver>
527  {
528  public:
537  ISTLBackend_OVLP_ILUn_Base (const GFS& gfs_, const C& c_, int n_=1, unsigned maxiter_=5000, int verbose_=1)
538  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), n(n_), maxiter(maxiter_), verbose(verbose_)
539  {}
540 
548  template<class M, class V, class W>
549  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
550  {
551  using Backend::Native;
552  using Backend::native;
553  typedef OverlappingOperator<C,M,V,W> POP;
554  POP pop(c,A);
555  typedef OVLPScalarProduct<GFS,V> PSP;
556  PSP psp(*this);
557  typedef SeqILUn<
558  Native<M>,
559  Native<V>,
560  Native<W>,
561  1
562  > SeqPrec;
563  SeqPrec seqprec(native(A),n,1.0);
565  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
566  int verb=0;
567  if (gfs.gridView().comm().rank()==0) verb=verbose;
568  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
569  Dune::InverseOperatorResult stat;
570  solver.apply(z,r,stat);
571  res.converged = stat.converged;
572  res.iterations = stat.iterations;
573  res.elapsed = stat.elapsed;
574  res.reduction = stat.reduction;
575  res.conv_rate = stat.conv_rate;
576  }
577  private:
578  const GFS& gfs;
579  const C& c;
580  int n;
581  unsigned maxiter;
582  int steps;
583  int verbose;
584  };
585 
588 
594  template<class GFS, class CC>
596  : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>
597  {
598  public:
607  ISTLBackend_OVLP_BCGS_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
608  int steps=5, int verbose=1)
609  : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>(gfs, cc, maxiter, steps, verbose)
610  {}
611  };
617  template<class GFS, class CC>
619  : public ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>
620  {
621  public:
629  ISTLBackend_OVLP_BCGS_ILU0 (const GFS& gfs, const CC& cc, unsigned maxiter=5000, int verbose=1)
630  : ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, maxiter, verbose)
631  {}
632  };
638  template<class GFS, class CC>
640  : public ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>
641  {
642  public:
651  ISTLBackend_OVLP_BCGS_ILUn (const GFS& gfs, const CC& cc, int n=1, unsigned maxiter=5000, int verbose=1)
652  : ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, n, maxiter, verbose)
653  {}
654  };
660  template<class GFS, class CC>
662  : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>
663  {
664  public:
673  ISTLBackend_OVLP_CG_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
674  int steps=5, int verbose=1)
675  : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>(gfs, cc, maxiter, steps, verbose)
676  {}
677  };
678 
684  template<class GFS, class CC>
687  {
688  public:
696  ISTLBackend_OVLP_GMRES_ILU0 (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000, int verbose_=1,
697  int restart_ = 20)
698  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), cc(cc_), maxiter(maxiter_), verbose(verbose_),
699  restart(restart_)
700  {}
701 
708  template<class M, class V, class W>
709  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
710  {
711  using Backend::Native;
712  using Backend::native;
713  typedef OverlappingOperator<CC,M,V,W> POP;
714  POP pop(cc,A);
715  typedef OVLPScalarProduct<GFS,V> PSP;
716  PSP psp(*this);
717  typedef SeqILU0<
718  Native<M>,
719  Native<V>,
720  Native<W>,
721  1
722  > SeqPrec;
723  SeqPrec seqprec(native(A),1.0);
725  WPREC wprec(gfs,seqprec,cc,this->parallelHelper());
726  int verb=0;
727  if (gfs.gridView().comm().rank()==0) verb=verbose;
728  RestartedGMResSolver<V> solver(pop,psp,wprec,reduction,restart,maxiter,verb);
729  Dune::InverseOperatorResult stat;
730  solver.apply(z,r,stat);
731  res.converged = stat.converged;
732  res.iterations = stat.iterations;
733  res.elapsed = stat.elapsed;
734  res.reduction = stat.reduction;
735  res.conv_rate = stat.conv_rate;
736  }
737 
738  private:
739  const GFS& gfs;
740  const CC& cc;
741  unsigned maxiter;
742  int steps;
743  int verbose;
744  int restart;
745  };
746 
748 
749  template<class GFS, class C, template<typename> class Solver>
752  {
753  public:
761  ISTLBackend_OVLP_SuperLU_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
762  int verbose_=1)
763  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
764  {}
765 
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  typedef OverlappingOperator<C,M,V,W> POP;
777  POP pop(c,A);
778  typedef OVLPScalarProduct<GFS,V> PSP;
779  PSP psp(*this);
780 #if HAVE_SUPERLU
781  typedef SuperLUSubdomainSolver<GFS,M,V,W> PREC;
782  PREC prec(gfs,A);
783  int verb=0;
784  if (gfs.gridView().comm().rank()==0) verb=verbose;
785  Solver<V> solver(pop,psp,prec,reduction,maxiter,verb);
786  Dune::InverseOperatorResult stat;
787  solver.apply(z,r,stat);
788  res.converged = stat.converged;
789  res.iterations = stat.iterations;
790  res.elapsed = stat.elapsed;
791  res.reduction = stat.reduction;
792  res.conv_rate = stat.conv_rate;
793 #else
794  std::cout << "No superLU support, please install and configure it." << std::endl;
795 #endif
796  }
797 
798  private:
799  const GFS& gfs;
800  const C& c;
801  unsigned maxiter;
802  int verbose;
803  };
804 
807 
812  template<class GFS, class CC>
814  : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>
815  {
816  public:
817 
825  ISTLBackend_OVLP_BCGS_SuperLU (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000,
826  int verbose_=1)
827  : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs_,cc_,maxiter_,verbose_)
828  {}
829  };
830 
836  template<class GFS, class CC>
838  : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>
839  {
840  public:
841 
849  ISTLBackend_OVLP_CG_SuperLU (const GFS& gfs_, const CC& cc_,
850  unsigned maxiter_=5000,
851  int verbose_=1)
852  : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
853  {}
854  };
855 
856 
860  template<class GFS>
862  : public LinearResultStorage
863  {
864  public:
869  explicit ISTLBackend_OVLP_ExplicitDiagonal (const GFS& gfs_)
870  : gfs(gfs_)
871  {}
872 
874  : gfs(other_.gfs)
875  {}
876 
881  template<class V>
882  typename V::ElementType norm(const V& v) const
883  {
884  static_assert
886  "ISTLBackend_OVLP_ExplicitDiagonal::norm() should not be "
887  "neccessary, so we skipped the implementation. If you have a "
888  "scenario where you need it, please implement it or report back to "
889  "us.");
890  }
891 
899  template<class M, class V, class W>
900  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
901  {
902  using Backend::Native;
903  using Backend::native;
904  Dune::SeqJac<
905  Native<M>,
906  Native<V>,
907  Native<W>
908  > jac(native(A),1,1.0);
909  jac.pre(native(z),native(r));
910  jac.apply(native(z),native(r));
911  jac.post(native(z));
912  if (gfs.gridView().comm().size()>1)
913  {
914  CopyDataHandle<GFS,V> copydh(gfs,z);
915  gfs.gridView().communicate(copydh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
916  }
917  res.converged = true;
918  res.iterations = 1;
919  res.elapsed = 0.0;
920  res.reduction = static_cast<double>(reduction);
921  res.conv_rate = static_cast<double>(reduction); // pow(reduction,1.0/1)
922  }
923 
924  private:
925  const GFS& gfs;
926  };
928 
929  template<class GO, int s, template<class,class,class,int> class Preconditioner,
930  template<class> class Solver>
932  {
933  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
935  typedef typename GO::Traits::Jacobian M;
936  typedef Backend::Native<M> MatrixType;
937  typedef typename GO::Traits::Domain V;
938  typedef Backend::Native<V> VectorType;
940 #if HAVE_MPI
941  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
942  typedef Dune::BlockPreconditioner<VectorType,VectorType,Comm,Smoother> ParSmoother;
943  typedef Dune::OverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
944 #else
945  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
946  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
947 #endif
948  typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
949  typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
950 
951  typedef typename V::ElementType RF;
952 
953  public:
954 
958  typedef Dune::Amg::Parameters Parameters;
959 
960  public:
961  ISTLBackend_AMG(const GFS& gfs_, unsigned maxiter_=5000,
962  int verbose_=1, bool reuse_=false,
963  bool usesuperlu_=true)
964  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), params(15,2000),
965  verbose(verbose_), reuse(reuse_), firstapply(true),
966  usesuperlu(usesuperlu_)
967  {
968  params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
969  params.setDebugLevel(verbose_);
970 #if !HAVE_SUPERLU
971  if (gfs.gridView().comm().rank() == 0 && usesuperlu == true)
972  {
973  std::cout << "WARNING: You are using AMG without SuperLU!"
974  << " Please consider installing SuperLU,"
975  << " or set the usesuperlu flag to false"
976  << " to suppress this warning." << std::endl;
977  }
978 #endif
979  }
980 
985  void setParameters(const Parameters& params_)
986  {
987  params = params_;
988  }
989 
990  void setparams(Parameters params_) DUNE_DEPRECATED_MSG("setparams() is deprecated, use setParameters() instead")
991  {
992  params = params_;
993  }
994 
1002  const Parameters& parameters() const
1003  {
1004  return params;
1005  }
1006 
1008  void setReuse(bool reuse_)
1009  {
1010  reuse = reuse_;
1011  }
1012 
1014  bool getReuse() const
1015  {
1016  return reuse;
1017  }
1018 
1023  typename V::ElementType norm (const V& v) const
1024  {
1025  typedef OverlappingScalarProduct<GFS,V> PSP;
1026  PSP psp(gfs,phelper);
1027  return psp.norm(v);
1028  }
1029 
1037  void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
1038  {
1039  Timer watch;
1040  Comm oocc(gfs.gridView().comm());
1041  MatrixType& mat=Backend::native(A);
1042  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
1043  Dune::Amg::FirstDiagonal> > Criterion;
1044 #if HAVE_MPI
1045  phelper.createIndexSetAndProjectForAMG(A, oocc);
1046  Operator oop(mat, oocc);
1047  Dune::OverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
1048 #else
1049  Operator oop(mat);
1050  Dune::SeqScalarProduct<VectorType> sp;
1051 #endif
1052  SmootherArgs smootherArgs;
1053  smootherArgs.iterations = 1;
1054  smootherArgs.relaxationFactor = 1;
1055  Criterion criterion(params);
1056  stats.tprepare=watch.elapsed();
1057  watch.reset();
1058 
1059  int verb=0;
1060  if (gfs.gridView().comm().rank()==0) verb=verbose;
1061  //only construct a new AMG if the matrix changes
1062  if (reuse==false || firstapply==true){
1063  amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
1064  firstapply = false;
1065  stats.tsetup = watch.elapsed();
1066  stats.levels = amg->maxlevels();
1067  stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
1068  }
1069  watch.reset();
1070  Solver<VectorType> solver(oop,sp,*amg,RF(reduction),maxiter,verb);
1071  Dune::InverseOperatorResult stat;
1072 
1073  solver.apply(Backend::native(z),Backend::native(r),stat);
1074  stats.tsolve= watch.elapsed();
1075  res.converged = stat.converged;
1076  res.iterations = stat.iterations;
1077  res.elapsed = stat.elapsed;
1078  res.reduction = stat.reduction;
1079  res.conv_rate = stat.conv_rate;
1080  }
1081 
1087  {
1088  return stats;
1089  }
1090 
1091  private:
1092  const GFS& gfs;
1093  PHELPER phelper;
1094  unsigned maxiter;
1095  Parameters params;
1096  int verbose;
1097  bool reuse;
1098  bool firstapply;
1099  bool usesuperlu;
1100  shared_ptr<AMG> amg;
1101  ISTLAMGStatistics stats;
1102  };
1103 
1106 
1113  template<class GO, int s=96>
1115  : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1116  {
1117  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1118  public:
1128  ISTLBackend_CG_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1129  int verbose_=1, bool reuse_=false,
1130  bool usesuperlu_=true)
1131  : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1132  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1133  {}
1134  };
1135 
1142  template<class GO, int s=96>
1144  : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1145  {
1146  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1147  public:
1157  ISTLBackend_BCGS_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1158  int verbose_=1, bool reuse_=false,
1159  bool usesuperlu_=true)
1160  : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1161  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1162  {}
1163  };
1164 
1171  template<class GO, int s=96>
1173  : public ISTLBackend_AMG<GO, s, Dune::SeqILU0, Dune::BiCGSTABSolver>
1174  {
1175  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1176  public:
1186  ISTLBackend_BCGS_AMG_ILU0(const GFS& gfs_, unsigned maxiter_=5000,
1187  int verbose_=1, bool reuse_=false,
1188  bool usesuperlu_=true)
1189  : ISTLBackend_AMG<GO, s, Dune::SeqILU0, Dune::BiCGSTABSolver>
1190  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1191  {}
1192  };
1193 
1195 
1197 
1198  } // namespace PDELab
1199 } // namespace Dune
1200 
1201 #endif
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:199
Definition: istl/ovlpistlsolverbackend.hh:323
STL namespace.
Definition: solver.hh:42
Overlapping parallel BiCGStab solver with SuperLU preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:813
virtual double norm(const X &x)
Norm of a right-hand side vector. The vector must be consistent on the interior+border partition...
Definition: istl/ovlpistlsolverbackend.hh:122
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: istl/ovlpistlsolverbackend.hh:607
OverlappingScalarProduct(const GFS &gfs_, const istl::ParallelHelper< GFS > &helper_)
Constructor needs to know the grid function space.
Definition: istl/ovlpistlsolverbackend.hh:101
ISTLBackend_BCGS_AMG_ILU0(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: istl/ovlpistlsolverbackend.hh:1186
X::ElementType field_type
Definition: istl/ovlpistlsolverbackend.hh:49
Definition: parallelhelper.hh:51
ISTLBackend_CG_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: istl/ovlpistlsolverbackend.hh:1128
ISTLBackend_OVLP_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int steps_=5, int verbose_=1)
make a linear solver object
Definition: istl/ovlpistlsolverbackend.hh:411
ISTLBackend_BCGS_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: istl/ovlpistlsolverbackend.hh:1157
void apply(M &A, V &z, V &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: istl/ovlpistlsolverbackend.hh:1037
ISTLBackend_OVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: istl/ovlpistlsolverbackend.hh:869
static const unsigned int value
Definition: gridfunctionspace/tags.hh:177
virtual void applyscaleadd(field_type alpha, const domain_type &x, range_type &y) const
apply operator to x, scale and add:
Definition: istl/ovlpistlsolverbackend.hh:67
Definition: istl/ovlpistlsolverbackend.hh:52
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper...
Definition: backend/interface.hh:183
Definition: istl/ovlpistlsolverbackend.hh:931
Dune::PDELab::Backend::Vector< GFS, typename P::domain_type::field_type > domain_type
The domain type of the preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:140
X::ElementType dot(const X &x, const X &y) const
Dot product of two vectors. It is assumed that the vectors are consistent on the interior+border part...
Definition: istl/ovlpistlsolverbackend.hh:335
Definition: genericdatahandle.hh:685
Definition: istl/ovlpistlsolverbackend.hh:399
OverlappingOperator(const CC &cc_, const M &A)
Definition: istl/ovlpistlsolverbackend.hh:54
OVLPScalarProductImplementation(const GFS &gfs_)
Definition: istl/ovlpistlsolverbackend.hh:326
Definition: genericdatahandle.hh:622
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: istl/ovlpistlsolverbackend.hh:774
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR.
Definition: istl/ovlpistlsolverbackend.hh:1143
virtual X::BaseT::field_type dot(const X &x, const X &y)
Definition: istl/ovlpistlsolverbackend.hh:381
ISTLBackend_OVLP_ExplicitDiagonal(const ISTLBackend_OVLP_ExplicitDiagonal &other_)
Definition: istl/ovlpistlsolverbackend.hh:873
Definition: istl/ovlpistlsolverbackend.hh:372
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: istl/ovlpistlsolverbackend.hh:709
ISTLBackend_OVLP_BCGS_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: istl/ovlpistlsolverbackend.hh:825
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: istl/ovlpistlsolverbackend.hh:1002
Overlapping parallel BiCGStab solver with SSOR preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:595
ISTLBackend_AMG(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: istl/ovlpistlsolverbackend.hh:961
ISTLBackend_OVLP_SuperLU_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: istl/ovlpistlsolverbackend.hh:761
OverlappingWrappedPreconditioner(const GFS &gfs_, P &prec_, const CC &cc_, const istl::ParallelHelper< GFS > &helper_)
Constructor.
Definition: istl/ovlpistlsolverbackend.hh:151
Definition: adaptivity.hh:27
Dune::PDELab::Backend::Vector< GFS, typename P::range_type::field_type > range_type
The range type of the preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:142
Definition: istl/ovlpistlsolverbackend.hh:88
Overlapping parallel CG solver with SuperLU preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:837
M matrix_type
export types
Definition: istl/ovlpistlsolverbackend.hh:46
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:415
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:798
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: istl/ovlpistlsolverbackend.hh:1008
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: istl/ovlpistlsolverbackend.hh:424
Definition: istl/ovlpistlsolverbackend.hh:750
void setparams(Parameters params_)
Definition: istl/ovlpistlsolverbackend.hh:990
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: istl/ovlpistlsolverbackend.hh:651
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: istl/ovlpistlsolverbackend.hh:549
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: istl/ovlpistlsolverbackend.hh:1014
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: istl/ovlpistlsolverbackend.hh:696
ISTLBackend_OVLP_ILUn_Base(const GFS &gfs_, const C &c_, int n_=1, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: istl/ovlpistlsolverbackend.hh:537
X domain_type
export types
Definition: istl/ovlpistlsolverbackend.hh:93
virtual void apply(const domain_type &x, range_type &y) const
apply operator to x:
Definition: istl/ovlpistlsolverbackend.hh:59
Dune::Amg::Parameters Parameters
Parameters object to customize matrix hierachy building.
Definition: istl/ovlpistlsolverbackend.hh:958
const istl::ParallelHelper< GFS > & parallelHelper() const
Definition: istl/ovlpistlsolverbackend.hh:354
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: istl/ovlpistlsolverbackend.hh:1023
X domain_type
Definition: istl/ovlpistlsolverbackend.hh:47
Class providing some statistics of the AMG solver.
Definition: istl/seqistlsolverbackend.hh:541
Definition: istl/ovlpistlsolverbackend.hh:41
virtual X::BaseT::field_type norm(const X &x)
Definition: istl/ovlpistlsolverbackend.hh:386
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by ILU0.
Definition: istl/ovlpistlsolverbackend.hh:1172
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:639
const ISTLAMGStatistics & statistics() const
Get statistics of the AMG solver (no of levels, timings).
Definition: istl/ovlpistlsolverbackend.hh:1086
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: istl/ovlpistlsolverbackend.hh:882
Dune::template FieldTraits< typename X::ElementType >::real_type norm(const X &x) const
Norm of a right-hand side vector. The vector must be consistent on the interior+border partition...
Definition: istl/ovlpistlsolverbackend.hh:348
virtual field_type dot(const X &x, const X &y)
Dot product of two vectors. It is assumed that the vectors are consistent on the interior+border part...
Definition: istl/ovlpistlsolverbackend.hh:110
Overlapping parallel CGS solver with SSOR preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:661
virtual const M & getmat() const
get matrix via *
Definition: istl/ovlpistlsolverbackend.hh:75
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:167
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:113
Definition: istl/ovlpistlsolverbackend.hh:525
ISTLBackend_OVLP_CG_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: istl/ovlpistlsolverbackend.hh:849
OVLPScalarProduct(const OVLPScalarProductImplementation< GFS > &implementation_)
Definition: istl/ovlpistlsolverbackend.hh:377
Overlapping parallel restarted GMRes solver with ILU0 preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:685
const std::string s
Definition: function.hh:1102
istl::ParallelHelper< GFS > & parallelHelper()
Definition: istl/ovlpistlsolverbackend.hh:360
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename W::ElementType >::real_type reduction)
solve the given linear system
Definition: istl/ovlpistlsolverbackend.hh:900
void setParameters(const Parameters &params_)
set AMG parameters
Definition: istl/ovlpistlsolverbackend.hh:985
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: istl/ovlpistlsolverbackend.hh:673
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:618
ISTLBackend_OVLP_BCGS_ILU0(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int verbose=1)
make a linear solver object
Definition: istl/ovlpistlsolverbackend.hh:629
Overlapping parallel conjugate gradient solver preconditioned with AMG smoothed by SSOR...
Definition: istl/ovlpistlsolverbackend.hh:1114
ISTLBackend_OVLP_ILU0_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: istl/ovlpistlsolverbackend.hh:474
virtual void post(domain_type &x)
Clean up.
Definition: istl/ovlpistlsolverbackend.hh:180
Y range_type
Definition: istl/ovlpistlsolverbackend.hh:48
X::ElementType field_type
Definition: istl/ovlpistlsolverbackend.hh:94
virtual void pre(domain_type &x, range_type &b)
Prepare the preconditioner.
Definition: istl/ovlpistlsolverbackend.hh:159
Definition: istl/ovlpistlsolverbackend.hh:463
Definition: istl/ovlpistlsolverbackend.hh:134
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: istl/ovlpistlsolverbackend.hh:486
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: istl/ovlpistlsolverbackend.hh:861