dune-pdelab  2.0.0
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  istl::raw(_A_).mv(istl::raw(x),istl::raw(y));
63  }
64 
66  virtual void applyscaleadd (field_type alpha, const domain_type& x, range_type& y) const
67  {
68  istl::raw(_A_).usmv(alpha,istl::raw(x),istl::raw(y));
70  }
71 
73  virtual const M& getmat () const
74  {
75  return _A_;
76  }
77 
78  private:
79  const CC& cc;
80  const M& _A_;
81  };
82 
83  // new scalar product assuming at least overlap 1
84  // uses unique partitioning of nodes for parallelization
85  template<class GFS, class X>
87  : public Dune::ScalarProduct<X>
88  {
89  public:
91  typedef X domain_type;
92  typedef typename X::ElementType field_type;
93 
95  enum {category=Dune::SolverCategory::overlapping};
96 
99  OverlappingScalarProduct (const GFS& gfs_, const istl::ParallelHelper<GFS>& helper_)
100  : gfs(gfs_), helper(helper_)
101  {}
102 
103 
108  virtual field_type dot (const X& x, const X& y)
109  {
110  // do local scalar product on unique partition
111  field_type sum = helper.disjointDot(x,y);
112 
113  // do global communication
114  return gfs.gridView().comm().sum(sum);
115  }
116 
120  virtual double norm (const X& x)
121  {
122  return sqrt(static_cast<double>(this->dot(x,x)));
123  }
124 
125  private:
126  const GFS& gfs;
127  const istl::ParallelHelper<GFS>& helper;
128  };
129 
130  // wrapped sequential preconditioner
131  template<class CC, class GFS, class P>
133  : public Dune::Preconditioner<typename Dune::PDELab::BackendVectorSelector<GFS,typename P::domain_type::field_type>::Type,
134  typename Dune::PDELab::BackendVectorSelector<GFS,typename P::range_type::field_type>::Type>
135  {
136  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(istl::raw(v),istl::raw(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(istl::raw(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 typename M::BaseT 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(istl::raw(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(istl::raw(v),istl::raw(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(istl::raw(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  Dune::InverseOperatorResult stat;
299  Y b(d); // need copy, since solver overwrites right hand side
300  solver.apply(istl::raw(v),istl::raw(b),stat);
301  if (gfs.gridView().comm().size()>1)
302  {
303  helper.maskForeignDOFs(istl::raw(v));
304  AddDataHandle<GFS,X> adddh(gfs,v);
305  gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
306  }
307  }
308 
312  virtual void post (X& x) {}
313 
314  private:
315  const GFS& gfs;
316  Dune::SuperLU<ISTLM> solver;
317  const istl::ParallelHelper<GFS>& helper;
318  };
319 #endif
320 
321  template<typename GFS>
323  {
324  public:
326  : gfs(gfs_), helper(gfs_)
327  {}
328 
333  template<typename X>
334  typename X::ElementType dot (const X& x, const X& y) const
335  {
336  // do local scalar product on unique partition
337  typename X::ElementType sum = helper.disjointDot(x,y);
338 
339  // do global communication
340  return gfs.gridView().comm().sum(sum);
341  }
342 
346  template<typename X>
347  typename X::ElementType norm (const X& x) const
348  {
349  return sqrt(static_cast<double>(this->dot(x,x)));
350  }
351 
353  {
354  return helper;
355  }
356 
357  // need also non-const version;
358  istl::ParallelHelper<GFS>& parallelHelper() // P.B.: needed for createIndexSetAndProjectForAMG
359  {
360  return helper;
361  }
362 
363  private:
364  const GFS& gfs;
366  };
367 
368 
369  template<typename GFS, typename X>
371  : public ScalarProduct<X>
372  {
373  public:
374  enum {category=Dune::SolverCategory::overlapping};
376  : implementation(implementation_)
377  {}
378 
379  virtual typename X::BaseT::field_type dot(const X& x, const X& y)
380  {
381  return implementation.dot(x,y);
382  }
383 
384  virtual typename X::BaseT::field_type norm (const X& x)
385  {
386  return sqrt(static_cast<double>(this->dot(x,x)));
387  }
388 
389  private:
390  const OVLPScalarProductImplementation<GFS>& implementation;
391  };
392 
393  template<class GFS, class C,
394  template<class,class,class,int> class Preconditioner,
395  template<class> class Solver>
398  {
399  public:
408  ISTLBackend_OVLP_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
409  int steps_=5, int verbose_=1)
410  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), steps(steps_), verbose(verbose_)
411  {}
412 
420  template<class M, class V, class W>
421  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
422  {
423  typedef OverlappingOperator<C,M,V,W> POP;
424  POP pop(c,A);
425  typedef OVLPScalarProduct<GFS,V> PSP;
426  PSP psp(*this);
427  typedef Preconditioner<typename M::BaseT,typename V::BaseT,typename W::BaseT,1> SeqPrec;
428  SeqPrec seqprec(istl::raw(A),steps,1.0);
430  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
431  int verb=0;
432  if (gfs.gridView().comm().rank()==0) verb=verbose;
433  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
434  Dune::InverseOperatorResult stat;
435  solver.apply(z,r,stat);
436  res.converged = stat.converged;
437  res.iterations = stat.iterations;
438  res.elapsed = stat.elapsed;
439  res.reduction = stat.reduction;
440  res.conv_rate = stat.conv_rate;
441  }
442  private:
443  const GFS& gfs;
444  const C& c;
445  unsigned maxiter;
446  int steps;
447  int verbose;
448  };
449 
450  // Base class for ILU0 as preconditioner
451  template<class GFS, class C,
452  template<class> class Solver>
455  {
456  public:
464  ISTLBackend_OVLP_ILU0_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000, int verbose_=1)
465  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
466  {}
467 
475  template<class M, class V, class W>
476  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
477  {
478  typedef OverlappingOperator<C,M,V,W> POP;
479  POP pop(c,A);
480  typedef OVLPScalarProduct<GFS,V> PSP;
481  PSP psp(*this);
482  typedef SeqILU0<typename M::BaseT,typename V::BaseT,typename W::BaseT,1> SeqPrec;
483  SeqPrec seqprec(istl::raw(A),1.0);
485  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
486  int verb=0;
487  if (gfs.gridView().comm().rank()==0) verb=verbose;
488  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
489  Dune::InverseOperatorResult stat;
490  solver.apply(z,r,stat);
491  res.converged = stat.converged;
492  res.iterations = stat.iterations;
493  res.elapsed = stat.elapsed;
494  res.reduction = stat.reduction;
495  res.conv_rate = stat.conv_rate;
496  }
497  private:
498  const GFS& gfs;
499  const C& c;
500  unsigned maxiter;
501  int steps;
502  int verbose;
503  };
504 
505  // Base class for ILUn as preconditioner
506  template<class GFS, class C,
507  template<class> class Solver>
510  {
511  public:
520  ISTLBackend_OVLP_ILUn_Base (const GFS& gfs_, const C& c_, int n_=1, unsigned maxiter_=5000, int verbose_=1)
521  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), n(n_), maxiter(maxiter_), verbose(verbose_)
522  {}
523 
531  template<class M, class V, class W>
532  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
533  {
534  typedef OverlappingOperator<C,M,V,W> POP;
535  POP pop(c,A);
536  typedef OVLPScalarProduct<GFS,V> PSP;
537  PSP psp(*this);
538  typedef SeqILUn<typename M::BaseT,typename V::BaseT,typename W::BaseT,1> SeqPrec;
539  SeqPrec seqprec(istl::raw(A),n,1.0);
541  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
542  int verb=0;
543  if (gfs.gridView().comm().rank()==0) verb=verbose;
544  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
545  Dune::InverseOperatorResult stat;
546  solver.apply(z,r,stat);
547  res.converged = stat.converged;
548  res.iterations = stat.iterations;
549  res.elapsed = stat.elapsed;
550  res.reduction = stat.reduction;
551  res.conv_rate = stat.conv_rate;
552  }
553  private:
554  const GFS& gfs;
555  const C& c;
556  int n;
557  unsigned maxiter;
558  int steps;
559  int verbose;
560  };
561 
564 
570  template<class GFS, class CC>
572  : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>
573  {
574  public:
583  ISTLBackend_OVLP_BCGS_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
584  int steps=5, int verbose=1)
585  : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>(gfs, cc, maxiter, steps, verbose)
586  {}
587  };
593  template<class GFS, class CC>
595  : public ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>
596  {
597  public:
605  ISTLBackend_OVLP_BCGS_ILU0 (const GFS& gfs, const CC& cc, unsigned maxiter=5000, int verbose=1)
606  : ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, maxiter, verbose)
607  {}
608  };
614  template<class GFS, class CC>
616  : public ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>
617  {
618  public:
627  ISTLBackend_OVLP_BCGS_ILUn (const GFS& gfs, const CC& cc, int n=1, unsigned maxiter=5000, int verbose=1)
628  : ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, n, maxiter, verbose)
629  {}
630  };
636  template<class GFS, class CC>
638  : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>
639  {
640  public:
649  ISTLBackend_OVLP_CG_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
650  int steps=5, int verbose=1)
651  : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>(gfs, cc, maxiter, steps, verbose)
652  {}
653  };
654 
660  template<class GFS, class CC>
663  {
664  public:
672  ISTLBackend_OVLP_GMRES_ILU0 (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000, int verbose_=1,
673  int restart_ = 20, bool recalc_defect_ = false)
674  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), cc(cc_), maxiter(maxiter_), verbose(verbose_),
675  restart(restart_), recalc_defect(recalc_defect_)
676  {}
677 
684  template<class M, class V, class W>
685  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
686  {
687  typedef OverlappingOperator<CC,M,V,W> POP;
688  POP pop(cc,A);
689  typedef OVLPScalarProduct<GFS,V> PSP;
690  PSP psp(*this);
691  typedef SeqILU0<typename M::BaseT,typename V::BaseT,typename W::BaseT,1> SeqPrec;
692  SeqPrec seqprec(istl::raw(A),1.0);
694  WPREC wprec(gfs,seqprec,cc,this->parallelHelper());
695  int verb=0;
696  if (gfs.gridView().comm().rank()==0) verb=verbose;
697  RestartedGMResSolver<V> solver(pop,psp,wprec,reduction,restart,maxiter,verb,recalc_defect);
698  Dune::InverseOperatorResult stat;
699  solver.apply(z,r,stat);
700  res.converged = stat.converged;
701  res.iterations = stat.iterations;
702  res.elapsed = stat.elapsed;
703  res.reduction = stat.reduction;
704  res.conv_rate = stat.conv_rate;
705  }
706 
707  private:
708  const GFS& gfs;
709  const CC& cc;
710  unsigned maxiter;
711  int steps;
712  int verbose;
713  int restart;
714  bool recalc_defect;
715  };
716 
718 
719  template<class GFS, class C, template<typename> class Solver>
722  {
723  public:
731  ISTLBackend_OVLP_SuperLU_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
732  int verbose_=1)
733  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
734  {}
735 
743  template<class M, class V, class W>
744  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
745  {
746  typedef OverlappingOperator<C,M,V,W> POP;
747  POP pop(c,A);
748  typedef OVLPScalarProduct<GFS,V> PSP;
749  PSP psp(*this);
750 #if HAVE_SUPERLU
751  typedef SuperLUSubdomainSolver<GFS,M,V,W> PREC;
752  PREC prec(gfs,A);
753  int verb=0;
754  if (gfs.gridView().comm().rank()==0) verb=verbose;
755  Solver<V> solver(pop,psp,prec,reduction,maxiter,verb);
756  Dune::InverseOperatorResult stat;
757  solver.apply(z,r,stat);
758  res.converged = stat.converged;
759  res.iterations = stat.iterations;
760  res.elapsed = stat.elapsed;
761  res.reduction = stat.reduction;
762  res.conv_rate = stat.conv_rate;
763 #else
764  std::cout << "No superLU support, please install and configure it." << std::endl;
765 #endif
766  }
767 
768  private:
769  const GFS& gfs;
770  const C& c;
771  unsigned maxiter;
772  int verbose;
773  };
774 
777 
782  template<class GFS, class CC>
784  : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>
785  {
786  public:
787 
795  ISTLBackend_OVLP_BCGS_SuperLU (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000,
796  int verbose_=1)
797  : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs_,cc_,maxiter_,verbose_)
798  {}
799  };
800 
806  template<class GFS, class CC>
808  : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>
809  {
810  public:
811 
819  ISTLBackend_OVLP_CG_SuperLU (const GFS& gfs_, const CC& cc_,
820  unsigned maxiter_=5000,
821  int verbose_=1)
822  : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
823  {}
824  };
825 
826 
830  template<class GFS>
832  : public LinearResultStorage
833  {
834  public:
839  explicit ISTLBackend_OVLP_ExplicitDiagonal (const GFS& gfs_)
840  : gfs(gfs_)
841  {}
842 
844  : gfs(other_.gfs)
845  {}
846 
851  template<class V>
852  typename V::ElementType norm(const V& v) const
853  {
854  dune_static_assert
856  "ISTLBackend_OVLP_ExplicitDiagonal::norm() should not be "
857  "neccessary, so we skipped the implementation. If you have a "
858  "scenario where you need it, please implement it or report back to "
859  "us.");
860  }
861 
869  template<class M, class V, class W>
870  void apply(M& A, V& z, W& r, typename W::ElementType reduction)
871  {
872  Dune::SeqJac<typename M::BaseT,typename V::BaseT,typename W::BaseT> jac(istl::raw(A),1,1.0);
873  jac.pre(istl::raw(z),istl::raw(r));
874  jac.apply(istl::raw(z),istl::raw(r));
875  jac.post(istl::raw(z));
876  if (gfs.gridView().comm().size()>1)
877  {
878  CopyDataHandle<GFS,V> copydh(gfs,z);
879  gfs.gridView().communicate(copydh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
880  }
881  res.converged = true;
882  res.iterations = 1;
883  res.elapsed = 0.0;
884  res.reduction = static_cast<double>(reduction);
885  res.conv_rate = static_cast<double>(reduction); // pow(reduction,1.0/1)
886  }
887 
888  private:
889  const GFS& gfs;
890  };
892 
893  template<class GO, int s, template<class,class,class,int> class Preconditioner,
894  template<class> class Solver>
896  {
897  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
899  typedef typename GO::Traits::Jacobian M;
900  typedef typename M::BaseT MatrixType;
901  typedef typename GO::Traits::Domain V;
902  typedef typename V::BaseT VectorType;
904 #if HAVE_MPI
905  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
906  typedef Dune::BlockPreconditioner<VectorType,VectorType,Comm,Smoother> ParSmoother;
907  typedef Dune::OverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
908 #else
909  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
910  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
911 #endif
912  typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
913  typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
914 
915  typedef typename V::ElementType RF;
916 
917  public:
918 
922  typedef Dune::Amg::Parameters Parameters;
923 
924  public:
925  ISTLBackend_AMG(const GFS& gfs_, unsigned maxiter_=5000,
926  int verbose_=1, bool reuse_=false,
927  bool usesuperlu_=true)
928  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), params(15,2000),
929  verbose(verbose_), reuse(reuse_), firstapply(true),
930  usesuperlu(usesuperlu_)
931  {
932  params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
933  params.setDebugLevel(verbose_);
934 #if !HAVE_SUPERLU
935  if (gfs.gridView().comm().rank() == 0 && usesuperlu == true)
936  {
937  std::cout << "WARNING: You are using AMG without SuperLU!"
938  << " Please consider installing SuperLU,"
939  << " or set the usesuperlu flag to false"
940  << " to suppress this warning." << std::endl;
941  }
942 #endif
943  }
944 
949  void setParameters(const Parameters& params_)
950  {
951  params = params_;
952  }
953 
954  void setparams(Parameters params_) DUNE_DEPRECATED_MSG("setparams() is deprecated, use setParameters() instead")
955  {
956  params = params_;
957  }
958 
966  const Parameters& parameters() const
967  {
968  return params;
969  }
970 
975  typename V::ElementType norm (const V& v) const
976  {
978  PSP psp(gfs,phelper);
979  return psp.norm(v);
980  }
981 
989  void apply(M& A, V& z, V& r, typename V::ElementType reduction)
990  {
991  Timer watch;
992  Comm oocc(gfs.gridView().comm());
993  MatrixType& mat=istl::raw(A);
994  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
995  Dune::Amg::FirstDiagonal> > Criterion;
996 #if HAVE_MPI
997  phelper.createIndexSetAndProjectForAMG(A, oocc);
998  Operator oop(mat, oocc);
999  Dune::OverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
1000 #else
1001  Operator oop(mat);
1002  Dune::SeqScalarProduct<VectorType> sp;
1003 #endif
1004  SmootherArgs smootherArgs;
1005  smootherArgs.iterations = 1;
1006  smootherArgs.relaxationFactor = 1;
1007  Criterion criterion(params);
1008  stats.tprepare=watch.elapsed();
1009  watch.reset();
1010 
1011  int verb=0;
1012  if (gfs.gridView().comm().rank()==0) verb=verbose;
1013  //only construct a new AMG if the matrix changes
1014  if (reuse==false || firstapply==true){
1015  amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
1016  firstapply = false;
1017  stats.tsetup = watch.elapsed();
1018  stats.levels = amg->maxlevels();
1019  stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
1020  }
1021  watch.reset();
1022  Solver<VectorType> solver(oop,sp,*amg,RF(reduction),maxiter,verb);
1023  Dune::InverseOperatorResult stat;
1024 
1025  solver.apply(istl::raw(z),istl::raw(r),stat);
1026  stats.tsolve= watch.elapsed();
1027  res.converged = stat.converged;
1028  res.iterations = stat.iterations;
1029  res.elapsed = stat.elapsed;
1030  res.reduction = stat.reduction;
1031  res.conv_rate = stat.conv_rate;
1032  }
1033 
1039  {
1040  return stats;
1041  }
1042 
1043  private:
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  shared_ptr<AMG> amg;
1053  ISTLAMGStatistics stats;
1054  };
1055 
1058 
1065  template<class GO, int s=96>
1067  : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1068  {
1069  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1070  public:
1080  ISTLBackend_CG_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1081  int verbose_=1, bool reuse_=false,
1082  bool usesuperlu_=true)
1083  : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1084  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1085  {}
1086  };
1087 
1094  template<class GO, int s=96>
1096  : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1097  {
1098  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1099  public:
1109  ISTLBackend_BCGS_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1110  int verbose_=1, bool reuse_=false,
1111  bool usesuperlu_=true)
1112  : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1113  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1114  {}
1115  };
1116 
1123  template<class GO, int s=96>
1125  : public ISTLBackend_AMG<GO, s, Dune::SeqILU0, Dune::BiCGSTABSolver>
1126  {
1127  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1128  public:
1138  ISTLBackend_BCGS_AMG_ILU0(const GFS& gfs_, unsigned maxiter_=5000,
1139  int verbose_=1, bool reuse_=false,
1140  bool usesuperlu_=true)
1141  : ISTLBackend_AMG<GO, s, Dune::SeqILU0, Dune::BiCGSTABSolver>
1142  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1143  {}
1144  };
1145 
1147 
1149 
1150  } // namespace PDELab
1151 } // namespace Dune
1152 
1153 #endif
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:421
X domain_type
Definition: ovlpistlsolverbackend.hh:47
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:685
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:532
ISTLBackend_BCGS_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1109
ISTLBackend_BCGS_AMG_ILU0(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1138
Y range_type
Definition: ovlpistlsolverbackend.hh:48
Definition: ovlpistlsolverbackend.hh:41
ISTLBackend_OVLP_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int steps_=5, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:408
Definition: ovlpistlsolverbackend.hh:720
Class providing some statistics of the AMG solver.
Definition: seqistlsolverbackend.hh:472
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:627
X domain_type
export types
Definition: ovlpistlsolverbackend.hh:91
Dune::PDELab::BackendVectorSelector< GFS, typename P::domain_type::field_type >::Type domain_type
The domain type of the preconditioner.
Definition: ovlpistlsolverbackend.hh:139
double tsetup
The time needed for building the AMG hierarchy (coarsening).
Definition: seqistlsolverbackend.hh:484
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: ovlpistlsolverbackend.hh:108
virtual void post(domain_type &x)
Clean up.
Definition: ovlpistlsolverbackend.hh:180
Definition: genericdatahandle.hh:686
void setParameters(const Parameters &params_)
set AMG parameters
Definition: ovlpistlsolverbackend.hh:949
bool converged
Definition: solver.hh:32
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by ILU0.
Definition: ovlpistlsolverbackend.hh:1124
OVLPScalarProduct(const OVLPScalarProductImplementation< GFS > &implementation_)
Definition: ovlpistlsolverbackend.hh:375
ISTLBackend_AMG(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: ovlpistlsolverbackend.hh:925
double tsolve
The time spent in solving the system (without building the hierarchy.
Definition: seqistlsolverbackend.hh:482
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:615
OverlappingWrappedPreconditioner(const GFS &gfs_, P &prec_, const CC &cc_, const istl::ParallelHelper< GFS > &helper_)
Constructor.
Definition: ovlpistlsolverbackend.hh:151
void setparams(Parameters params_)
Definition: ovlpistlsolverbackend.hh:954
double elapsed
Definition: solver.hh:34
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: ovlpistlsolverbackend.hh:334
The category the preconditioner is part of.
Definition: ovlpistlsolverbackend.hh:147
void apply(M &A, V &z, W &r, typename W::ElementType reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:870
Overlapping parallel CGS solver with SSOR preconditioner.
Definition: ovlpistlsolverbackend.hh:637
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: common/constraints.hh:772
int levels
the number of levels in the AMG hierarchy.
Definition: seqistlsolverbackend.hh:480
virtual const M & getmat() const
get matrix via *
Definition: ovlpistlsolverbackend.hh:73
virtual void applyscaleadd(field_type alpha, const domain_type &x, range_type &y) const
apply operator to x, scale and add:
Definition: ovlpistlsolverbackend.hh:66
Definition: ovlpistlsolverbackend.hh:508
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:476
ISTLBackend_OVLP_ILU0_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:464
X::ElementType field_type
Definition: ovlpistlsolverbackend.hh:92
Definition: ovlpistlsolverbackend.hh:95
OVLPScalarProductImplementation(const GFS &gfs_)
Definition: ovlpistlsolverbackend.hh:325
istl::ParallelHelper< GFS > & parallelHelper()
Definition: ovlpistlsolverbackend.hh:358
V & raw(V &v)
Returns the raw ISTL object associated with v, or v itself it is already an ISTL object.
Definition: backend/istl/utility.hh:26
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:52
const istl::ParallelHelper< GFS > & parallelHelper() const
Definition: ovlpistlsolverbackend.hh:352
Definition: ovlpistlsolverbackend.hh:374
Definition: ovlpistlsolverbackend.hh:132
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: ovlpistlsolverbackend.hh:966
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:594
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: ovlpistlsolverbackend.hh:852
ISTLBackend_OVLP_GMRES_ILU0(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1, int restart_=20, bool recalc_defect_=false)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:672
Definition: parallelhelper.hh:45
ISTLBackend_OVLP_BCGS_ILU0(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:605
Overlapping parallel conjugate gradient solver preconditioned with AMG smoothed by SSOR...
Definition: ovlpistlsolverbackend.hh:1066
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: ovlpistlsolverbackend.hh:975
void apply(M &A, V &z, V &r, typename V::ElementType reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:989
Definition: genericdatahandle.hh:623
virtual double norm(const X &x)
Norm of a right-hand side vector. The vector must be consistent on the interior+border partition...
Definition: ovlpistlsolverbackend.hh:120
ISTLBackend_OVLP_SuperLU_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:731
Definition: ovlpistlsolverbackend.hh:52
void createIndexSetAndProjectForAMG(MatrixType &m, Comm &c)
Makes the matrix consistent and creates the parallel information for AMG.
Overlapping parallel BiCGStab solver with SuperLU preconditioner.
Definition: ovlpistlsolverbackend.hh:783
RFType reduction
Definition: solver.hh:35
ISTLBackend_OVLP_BCGS_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:795
ISTLBackend_OVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:839
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
Definition: ovlpistlsolverbackend.hh:453
BackendVectorSelectorHelper< Backend, GridFunctionSpace, FieldType >::Type Type
Definition: backendselector.hh:14
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: ovlpistlsolverbackend.hh:831
virtual X::BaseT::field_type dot(const X &x, const X &y)
Definition: ovlpistlsolverbackend.hh:379
virtual void pre(domain_type &x, range_type &b)
Prepare the preconditioner.
Definition: ovlpistlsolverbackend.hh:159
Definition: solver.hh:42
virtual void apply(const domain_type &x, range_type &y) const
apply operator to x:
Definition: ovlpistlsolverbackend.hh:59
const ISTLAMGStatistics & statistics() const
Get statistics of the AMG solver (no of levels, timings).
Definition: ovlpistlsolverbackend.hh:1038
Definition: ovlpistlsolverbackend.hh:322
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: ovlpistlsolverbackend.hh:520
bool directCoarseLevelSolver
True if a direct solver was used on the coarset level.
Definition: seqistlsolverbackend.hh:488
Definition: ovlpistlsolverbackend.hh:396
ISTLBackend_OVLP_ExplicitDiagonal(const ISTLBackend_OVLP_ExplicitDiagonal &other_)
Definition: ovlpistlsolverbackend.hh:843
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:583
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition: ovlpistlsolverbackend.hh:167
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR.
Definition: ovlpistlsolverbackend.hh:1095
M matrix_type
export types
Definition: ovlpistlsolverbackend.hh:46
Dune::PDELab::BackendVectorSelector< GFS, typename P::range_type::field_type >::Type range_type
The range type of the preconditioner.
Definition: ovlpistlsolverbackend.hh:142
ISTLBackend_OVLP_CG_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:819
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:417
Dune::Amg::Parameters Parameters
Parameters object to customize matrix hierachy building.
Definition: ovlpistlsolverbackend.hh:922
X::ElementType norm(const X &x) const
Norm of a right-hand side vector. The vector must be consistent on the interior+border partition...
Definition: ovlpistlsolverbackend.hh:347
X::ElementType field_type
Definition: ovlpistlsolverbackend.hh:49
Definition: ovlpistlsolverbackend.hh:86
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:744
Overlapping parallel restarted GMRes solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:661
Definition: ovlpistlsolverbackend.hh:895
double tprepare
The needed for computing the parallel information and for adapting the linear system.
Definition: seqistlsolverbackend.hh:478
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:649
Overlapping parallel CG solver with SuperLU preconditioner.
Definition: ovlpistlsolverbackend.hh:807
OverlappingOperator(const CC &cc_, const M &A)
Definition: ovlpistlsolverbackend.hh:54
virtual X::BaseT::field_type norm(const X &x)
Definition: ovlpistlsolverbackend.hh:384
RFType conv_rate
Definition: solver.hh:36
Definition: ovlpistlsolverbackend.hh:370
OverlappingScalarProduct(const GFS &gfs_, const istl::ParallelHelper< GFS > &helper_)
Constructor needs to know the grid function space.
Definition: ovlpistlsolverbackend.hh:99
Overlapping parallel BiCGStab solver with SSOR preconditioner.
Definition: ovlpistlsolverbackend.hh:571
ISTLBackend_CG_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1080
const std::string s
Definition: function.hh:1103
unsigned int iterations
Definition: solver.hh:33