dune-pdelab  2.5-dev
ovlp_amg_dg_backend.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH
2 #define DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH
3 
4 // this is here for backwards compatibility and deprecation warnings, remove after 2.5.0
5 #include "ensureistlinclude.hh"
6 
7 #include <dune/common/parametertree.hh>
8 #include <dune/common/power.hh>
9 
10 #include <dune/istl/matrixmatrix.hh>
11 
12 #include <dune/grid/common/datahandleif.hh>
22 
23 namespace Dune {
24  namespace PDELab {
25  //***********************************************************
26  // A data handle / function to communicate matrix entries
27  // in the overlap. It turned out that it is actually not
28  // required to do that, but anyway it is there now.
29  //***********************************************************
30 
33  template<class GFS>
35  : public Dune::CommDataHandleIF<LocalGlobalMapDataHandle<GFS>,int>
36  {
37  public:
38  // some types from the grid view
39  typedef typename GFS::Traits::GridView GV;
40  typedef typename GV::IndexSet IndexSet;
41  typedef typename IndexSet::IndexType IndexType;
42  typedef typename GV::Grid Grid;
43  typedef typename Grid::Traits::GlobalIdSet GlobalIdSet;
44  typedef typename GlobalIdSet::IdType IdType;
45 
47  typedef int DataType;
48 
49  // map types
50  typedef std::map<IndexType,IdType> LocalToGlobalMap;
51  typedef std::map<IdType,IndexType> GlobalToLocalMap;
52 
54  bool contains (int dim, int codim) const
55  {
56  return (codim==0);
57  }
58 
60  bool fixedsize (int dim, int codim) const
61  {
62  return true;
63  }
64 
69  template<class EntityType>
70  size_t size (EntityType& e) const
71  {
72  return 1;
73  }
74 
76  template<class MessageBuffer, class EntityType>
77  void gather (MessageBuffer& buff, const EntityType& e) const
78  {
79  // fill map
80  IndexType myindex = indexSet.index(e);
81  IdType myid = globalIdSet.id(e);
82  l2g[myindex] = myid;
83  g2l[myid] = myindex;
84  //std::cout << gv.comm().rank() << ": gather local=" << myindex << " global=" << myid << std::endl;
85 
86  buff.write(0); // this is a dummy, we are not interested in the data
87  }
88 
93  template<class MessageBuffer, class EntityType>
94  void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
95  {
96  DataType x;
97  buff.read(x); // this is a dummy, we are not interested in the data
98 
99  IndexType myindex = indexSet.index(e);
100  IdType myid = globalIdSet.id(e);
101  l2g[myindex] = myid;
102  g2l[myid] = myindex;
103  //std::cout << gv.comm().rank() << ": scatter local=" << myindex << " global=" << myid << std::endl;
104  }
105 
107  LocalGlobalMapDataHandle (const GFS& gfs_, LocalToGlobalMap& l2g_, GlobalToLocalMap& g2l_)
108  : gfs(gfs_), gv(gfs.gridView()), indexSet(gv.indexSet()),
109  grid(gv.grid()), globalIdSet(grid.globalIdSet()),
110  l2g(l2g_), g2l(g2l_)
111  {
112  }
113 
114  private:
115  const GFS& gfs;
116  GV gv;
117  const IndexSet& indexSet;
118  const Grid& grid;
119  const GlobalIdSet& globalIdSet;
120  LocalToGlobalMap& l2g;
121  GlobalToLocalMap& g2l;
122  };
123 
124 
125  // A DataHandle class to exchange rows of a sparse matrix
126  template<class GFS, class M> // mapper type and vector type
128  : public Dune::CommDataHandleIF<MatrixExchangeDataHandle<GFS,M>,
129  std::pair<typename GFS::Traits::GridView::Grid::Traits::GlobalIdSet::IdType,
130  typename M::block_type> >
131  {
132  public:
133  // some types from the grid view
134  typedef typename GFS::Traits::GridView GV;
135  typedef typename GV::IndexSet IndexSet;
136  typedef typename IndexSet::IndexType IndexType;
137  typedef typename GV::Grid Grid;
138  typedef typename Grid::Traits::GlobalIdSet GlobalIdSet;
139  typedef typename GlobalIdSet::IdType IdType;
140 
141  // some types from the matrix
142  typedef typename M::block_type B;
143 
145  typedef std::pair<IdType,B> DataType;
146 
147  // map types
148  typedef std::map<IndexType,IdType> LocalToGlobalMap;
149  typedef std::map<IdType,IndexType> GlobalToLocalMap;
150 
152  bool contains (int dim, int codim) const
153  {
154  return (codim==0);
155  }
156 
158  bool fixedsize (int dim, int codim) const
159  {
160  return false;
161  }
162 
167  template<class EntityType>
168  size_t size (EntityType& e) const
169  {
170  IndexType myindex = indexSet.index(e);
171  typename M::row_type myrow = m[myindex];
172  typename M::row_type::iterator endit=myrow.end();
173  size_t count=0;
174  for (typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it)
175  {
176  typename LocalToGlobalMap::const_iterator find=l2g.find(it.index());
177  if (find!=l2g.end())
178  {
179  count++;
180  }
181  }
182  //std::cout << gv.comm().rank() << ": row=" << myindex << " size=" << count << std::endl;
183  return count;
184  }
185 
187  template<class MessageBuffer, class EntityType>
188  void gather (MessageBuffer& buff, const EntityType& e) const
189  {
190  IndexType myindex = indexSet.index(e);
191  //std::cout << gv.comm().rank() << ": begin gather local=" << myindex << std::endl;
192  typename M::row_type myrow = m[myindex];
193  typename M::row_type::iterator endit=myrow.end();
194  size_t count = 0;
195  for (typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it)
196  {
197  typename LocalToGlobalMap::const_iterator find=l2g.find(it.index());
198  if (find!=l2g.end())
199  {
200  buff.write(std::make_pair(find->second,*it));
201  //std::cout << gv.comm().rank() << ": ==> local=" << find->first << " global=" << find->second << std::endl;
202  count++;
203  }
204  }
205  //std::cout << gv.comm().rank() << ": end gather row=" << myindex << " size=" << count << std::endl;
206  }
207 
212  template<class MessageBuffer, class EntityType>
213  void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
214  {
215  IndexType myindex = indexSet.index(e);
216  std::cout << gv.comm().rank() << ": begin scatter local=" << myindex << " size=" << n << std::endl;
217  DataType x;
218  size_t count = 0;
219  for (size_t i=0; i<n; ++i)
220  {
221  buff.read(x);
222  std::cout << gv.comm().rank() << ": --> received global=" << x.first << std::endl;
223  typename GlobalToLocalMap::const_iterator find=g2l.find(x.first);
224  if (find!=g2l.end())
225  {
226  IndexType nbindex=find->second;
227  if (m.exists(myindex,nbindex))
228  {
229  m[myindex][nbindex] = x.second;
230  B block(x.second);
231  block -= m2[myindex][nbindex];
232  std::cout << gv.comm().rank() << ": compare i=" << myindex << " j=" << nbindex
233  << " norm=" << block.infinity_norm() << std::endl;
234 
235  count++;
236  //std::cout << gv.comm().rank() << ": --> local=" << find->first << " global=" << find->second << std::endl;
237  }
238  }
239  }
240  //std::cout << gv.comm().rank() << ": end scatter row=" << myindex << " size=" << count << std::endl;
241  }
242 
244  MatrixExchangeDataHandle (const GFS& gfs_, M& m_, const LocalToGlobalMap& l2g_, const GlobalToLocalMap& g2l_,M& m2_)
245  : gfs(gfs_), m(m_), gv(gfs.gridView()), indexSet(gv.indexSet()),
246  grid(gv.grid()), globalIdSet(grid.globalIdSet()),
247  l2g(l2g_), g2l(g2l_), m2(m2_)
248  {
249  }
250 
251  private:
252  const GFS& gfs;
253  M& m;
254  GV gv;
255  const IndexSet& indexSet;
256  const Grid& grid;
257  const GlobalIdSet& globalIdSet;
258  const LocalToGlobalMap& l2g;
259  const GlobalToLocalMap& g2l;
260  M& m2;
261  };
262 
263 
266  template<class GFS, class T, class A, int n, int m>
267  void restore_overlap_entries (const GFS& gfs, Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix,
268  Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix2)
269  {
270  typedef Dune::FieldMatrix<T,n,m> B;
271  typedef Dune::BCRSMatrix<B,A> M; // m is of type M
274 
275  // build up two maps to associate local indices and global ids
276  LocalToGlobalMap l2g;
277  GlobalToLocalMap g2l;
278  LocalGlobalMapDataHandle<GFS> lgdh(gfs,l2g,g2l);
279  if (gfs.gridView().comm().size()>1)
280  gfs.gridView().communicate(lgdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
281 
282  // exchange matrix entries
283  MatrixExchangeDataHandle<GFS,M> mexdh(gfs,matrix,l2g,g2l,matrix2);
284  if (gfs.gridView().comm().size()>1)
285  gfs.gridView().communicate(mexdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
286  }
287 
288  //***********************************************************
289  // The DG/AMG preconditioner in the overlapping case
290  //***********************************************************
291 
301  template<class DGGFS, class DGMatrix, class DGPrec, class DGCC,
302  class CGGFS, class CGPrec, class CGCC,
303  class P, class DGHelper, class Comm>
305  : public Dune::Preconditioner<Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::domain_type::field_type>,
306  Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::range_type::field_type>>
307  {
308  const DGGFS& dggfs;
309  DGMatrix& dgmatrix;
310  DGPrec& dgprec;
311  const DGCC& dgcc;
312  const CGGFS& cggfs;
313  CGPrec& cgprec;
314  const CGCC& cgcc;
315  P& p;
316  const DGHelper& dghelper;
317  const Comm& comm;
318  int n1,n2;
319 
320  public:
327 
328  // define the category
329  enum {
331  category=Dune::SolverCategory::overlapping
332  };
333 
341  OvlpDGAMGPrec (const DGGFS& dggfs_, DGMatrix& dgmatrix_, DGPrec& dgprec_, const DGCC& dgcc_,
342  const CGGFS& cggfs_, CGPrec& cgprec_, const CGCC& cgcc_, P& p_,
343  const DGHelper& dghelper_, const Comm& comm_, int n1_, int n2_)
344  : dggfs(dggfs_), dgmatrix(dgmatrix_), dgprec(dgprec_), dgcc(dgcc_),
345  cggfs(cggfs_), cgprec(cgprec_), cgcc(cgcc_), p(p_), dghelper(dghelper_),
346  comm(comm_), n1(n1_), n2(n2_)
347  {
348  }
349 
355  virtual void pre (V& x, W& b)
356  {
357  using Backend::native;
358  dgprec.pre(native(x),native(b));
359  CGW cgd(cggfs,0.0);
360  CGV cgv(cggfs,0.0);
361  cgprec.pre(native(cgv),native(cgd));
362  }
363 
369  virtual void apply (V& x, const W& b)
370  {
371  using Backend::native;
372  // need local copies to store defect and solution
373  W d(b);
375  V v(x); // only to get correct size
376 
377  // pre-smoothing on DG matrix
378  for (int i=0; i<n1; i++)
379  {
380  using Backend::native;
381  v = 0.0;
382  dgprec.apply(native(v),native(d));
384  if (dggfs.gridView().comm().size()>1)
385  dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
386  dgmatrix.mmv(native(v),native(d));
388  x += v;
389  }
390 
391  // restrict defect to CG subspace
392  dghelper.maskForeignDOFs(d); // DG defect is additive for overlap 1, but in case we use more
393  CGW cgd(cggfs,0.0);
394  p.mtv(native(d),native(cgd));
396  if (cggfs.gridView().comm().size()>1)
397  cggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication); // now we have consistent defect on coarse grid
399  comm.project(native(cgd));
400  CGV cgv(cggfs,0.0);
401 
402 
403  // call preconditioner
404  cgprec.apply(native(cgv),native(cgd));
405 
406  // prolongate correction
407  p.mv(native(cgv),native(v));
408  dgmatrix.mmv(native(v),native(d));
410  x += v;
411 
412  // post-smoothing on DG matrix
413  for (int i=0; i<n2; i++)
414  {
415  v = 0.0;
416  dgprec.apply(native(v),native(d));
418  if (dggfs.gridView().comm().size()>1)
419  dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
420  dgmatrix.mmv(native(v),native(d));
422  x += v;
423  }
424  }
425 
431  virtual void post (V& x)
432  {
433  dgprec.post(Backend::native(x));
434  CGV cgv(cggfs,0.0);
435  cgprec.post(Backend::native(cgv));
436  }
437  };
438 
439 //***********************************************************
440 // The DG/AMG linear solver backend in the overlapping case
441 //***********************************************************
442 
455 template<class DGGO, class DGCC, class CGGFS, class CGCC, class TransferLOP,
456  template<class,class,class,int> class DGPrec, template<class> class Solver, int s=96>
458  public Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>,
460 {
461 public:
462  // DG grid function space
463  typedef typename DGGO::Traits::TrialGridFunctionSpace GFS;
464 
465  // vectors and matrices on DG level
466  typedef typename DGGO::Traits::Jacobian M; // wrapped istl DG matrix
467  typedef typename DGGO::Traits::Domain V; // wrapped istl DG vector
468  typedef Backend::Native<M> Matrix; // istl DG matrix
469  typedef Backend::Native<V> Vector; // istl DG vector
470  typedef typename Vector::field_type field_type;
471 
472  // vectors and matrices on CG level
473  using CGV = Dune::PDELab::Backend::Vector<CGGFS,field_type>; // wrapped istl CG vector
474  typedef Backend::Native<CGV> CGVector; // istl CG vector
475 
476  // prolongation matrix
479  typedef TransferLOP CGTODGLOP; // local operator
481  typedef typename PGO::Jacobian PMatrix; // wrapped ISTL prolongation matrix
482  typedef Backend::Native<PMatrix> P; // ISTL prolongation matrix
483 
484  // CG subspace matrix
485  typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG;
486  typedef typename Dune::MatMultMatResult<PTADG,P>::type CGMatrix; // istl coarse space matrix
487 
488  // AMG in CG-subspace
490  typedef Dune::OverlappingSchwarzOperator<CGMatrix,CGVector,CGVector,Comm> ParCGOperator;
491  typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother;
492  typedef Dune::BlockPreconditioner<CGVector,CGVector,Comm,Smoother> ParSmoother;
493  typedef Dune::Amg::AMG<ParCGOperator,CGVector,ParSmoother,Comm> AMG;
494  typedef Dune::Amg::Parameters Parameters;
495 
496 private:
497 
498  const GFS& gfs;
499  DGGO& dggo;
500  const DGCC& dgcc;
501  CGGFS& cggfs;
502  const CGCC& cgcc;
503  std::shared_ptr<AMG> amg;
504  Parameters amg_parameters;
505  unsigned maxiter;
506  int verbose;
507  bool reuse;
508  bool firstapply;
509  bool usesuperlu;
510  std::size_t low_order_space_entries_per_row;
511 
512  CGTODGLOP cgtodglop; // local operator to assemble prolongation matrix
513  PGO pgo; // grid operator to assemble prolongation matrix
514  PMatrix pmatrix; // wrapped prolongation matrix
515 
518  class EmptyLop : public Dune::PDELab::NumericalJacobianApplyVolume<EmptyLop>,
522  {
523  };
524 
525  // grid operator with empty local operator for constraints assembly in parallel
527  // CG-subspace matrix
528  typedef typename CGGO::Jacobian CGM;
529  CGM acg;
530 
531 public:
532 
533  // access to prolongation matrix
535  {
536  return pmatrix;
537  }
538 
543  void setParameters(const Parameters& amg_parameters_)
544  {
545  amg_parameters = amg_parameters_;
546  }
547 
555  const Parameters& parameters() const
556  {
557  return amg_parameters;
558  }
559 
561  void setReuse(bool reuse_)
562  {
563  reuse = reuse_;
564  }
565 
567  bool getReuse() const
568  {
569  return reuse;
570  }
571 
574  ISTLBackend_OVLP_AMG_4_DG(DGGO& dggo_, const DGCC& dgcc_, CGGFS& cggfs_, const CGCC& cgcc_,
575  unsigned maxiter_=5000, int verbose_=1, bool reuse_=false,
576  bool usesuperlu_=true)
577  : Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>(dggo_.trialGridFunctionSpace())
578  , gfs(dggo_.trialGridFunctionSpace())
579  , dggo(dggo_)
580  , dgcc(dgcc_)
581  , cggfs(cggfs_)
582  , cgcc(cgcc_)
583  , amg_parameters(15,2000)
584  , maxiter(maxiter_)
585  , verbose(verbose_)
586  , reuse(reuse_)
587  , firstapply(true)
588  , usesuperlu(usesuperlu_)
589  , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
590  , cgtodglop()
591  , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
592  , pmatrix(pgo)
593  , acg(Backend::attached_container())
594  {
595  amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
596  amg_parameters.setDebugLevel(verbose_);
597 #if !HAVE_SUPERLU
598  if (usesuperlu == true)
599  {
600  if (gfs.gridView().comm().rank()==0)
601  std::cout << "WARNING: You are using AMG without SuperLU!"
602  << " Please consider installing SuperLU,"
603  << " or set the usesuperlu flag to false"
604  << " to suppress this warning." << std::endl;
605  }
606 #endif
607 
608  // assemble prolongation matrix; this will not change from one apply to the next
609  pmatrix = 0.0;
610  if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
611  CGV cgx(cggfs,0.0); // need vector to call jacobian
612  pgo.jacobian(cgx,pmatrix);
613  }
614 
615 
618  ISTLBackend_OVLP_AMG_4_DG(DGGO& dggo_, const DGCC& dgcc_, CGGFS& cggfs_, const CGCC& cgcc_,
619  const ParameterTree& params)
620  : Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>(dggo_.trialGridFunctionSpace())
621  , gfs(dggo_.trialGridFunctionSpace())
622  , dggo(dggo_)
623  , dgcc(dgcc_)
624  , cggfs(cggfs_)
625  , cgcc(cgcc_)
626  , maxiter(params.get<int>("max_iterations",5000))
627  , amg_parameters(15,2000)
628  , verbose(params.get<int>("verbose",1))
629  , reuse(params.get<bool>("reuse",false))
630  , firstapply(true)
631  , usesuperlu(params.get<bool>("use_superlu",true))
632  , low_order_space_entries_per_row(params.get<std::size_t>("low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power))
633  , cgtodglop()
634  , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
635  , pmatrix(pgo)
636  , acg(Backend::attached_container())
637  {
638  amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
639  amg_parameters.setDebugLevel(params.get<int>("verbose",1));
640 #if !HAVE_SUPERLU
641  if (usesuperlu == true)
642  {
643  if (gfs.gridView().comm().rank()==0)
644  std::cout << "WARNING: You are using AMG without SuperLU!"
645  << " Please consider installing SuperLU,"
646  << " or set the usesuperlu flag to false"
647  << " to suppress this warning." << std::endl;
648  }
649 #endif
650 
651  // assemble prolongation matrix; this will not change from one apply to the next
652  pmatrix = 0.0;
653  if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
654  CGV cgx(cggfs,0.0); // need vector to call jacobian
655  pgo.jacobian(cgx,pmatrix);
656  }
657 
658 
666  void apply (M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
667  {
668  using Backend::native;
669  // make operator and scalar product for overlapping solver
671  POP pop(dgcc,A);
673  PSP psp(*this);
674 
675  // compute CG matrix
676  // make grid operator with empty local operator => matrix data type and constraints assembly
677  EmptyLop emptylop;
678  CGGO cggo(cggfs,cgcc,cggfs,cgcc,emptylop,MBE(low_order_space_entries_per_row));
679 
680  // do triple matrix product ACG = P^T ADG P; this is purely local
681  Dune::Timer watch;
682  watch.reset();
683  // only do triple matrix product if the matrix changes
684  double triple_product_time = 0.0;
685  // no need to set acg here back to zero, this is done in matMultmat
686  if(reuse == false || firstapply == true) {
687  PTADG ptadg;
688  Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A)); // 1a
689  //Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A2)); // 1b
690  Dune::matMultMat(native(acg),ptadg,native(pmatrix));
691  triple_product_time = watch.elapsed();
692  if (verbose>0 && gfs.gridView().comm().rank()==0)
693  std::cout << "=== triple matrix product " << triple_product_time << " s" << std::endl;
694  //Dune::printmatrix(std::cout,native(acg),"triple product matrix","row",10,2);
695  CGV cgx(cggfs,0.0); // need vector to call jacobian
696  cggo.jacobian(cgx,acg); // insert trivial rows at processor boundaries
697  //std::cout << "CG constraints: " << cgcc.size() << " out of " << cggfs.globalSize() << std::endl;
698  }
699  else if(verbose>0 && gfs.gridView().comm().rank()==0)
700  std::cout << "=== reuse CG matrix, SKIPPING triple matrix product " << std::endl;
701 
702  // NOW we need to insert the processor boundary conditions in DG matrix
704  DGGOEmpty dggoempty(gfs,dgcc,gfs,dgcc,emptylop,MBE(1 << GFS::Traits::GridView::dimension));
705  dggoempty.jacobian(z,A);
706 
707  // and in the residual
709 
710  // now set up parallel AMG solver for the CG subspace
711  Comm oocc(gfs.gridView().comm());
713  CGHELPER cghelper(cggfs,2);
714  cghelper.createIndexSetAndProjectForAMG(acg,oocc);
715  ParCGOperator paroop(native(acg),oocc);
716  Dune::OverlappingSchwarzScalarProduct<CGVector,Comm> sp(oocc);
717 
718  typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
719  SmootherArgs smootherArgs;
720  smootherArgs.iterations = 1;
721  smootherArgs.relaxationFactor = 1.0;
722  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
723  Criterion criterion(amg_parameters);
724  watch.reset();
725 
726  // only construct a new AMG for the CG-subspace if the matrix changes
727  double amg_setup_time = 0.0;
728  if(reuse == false || firstapply == true) {
729  amg.reset(new AMG(paroop,criterion,smootherArgs,oocc));
730  firstapply = false;
731  amg_setup_time = watch.elapsed();
732  if (verbose>0 && gfs.gridView().comm().rank()==0)
733  std::cout << "=== AMG setup " <<amg_setup_time << " s" << std::endl;
734  }
735  else if (verbose>0 && gfs.gridView().comm().rank()==0)
736  std::cout << "=== reuse CG matrix, SKIPPING AMG setup " << std::endl;
737 
738  // set up hybrid DG/CG preconditioner
739  typedef DGPrec<Matrix,Vector,Vector,1> DGPrecType;
740  DGPrecType dgprec(native(A),1,0.92);
741  //DGPrecType dgprec(native(A),0.92);
744  HybridPrec hybridprec(gfs,native(A),dgprec,dgcc,cggfs,*amg,cgcc,native(pmatrix),
745  this->parallelHelper(),oocc,3,3);
746 
747  // /********/
748  // /* Test */
749  // /********/
750  // Solver<CGVector> testsolver(paroop,sp,amg,1e-8,100,2);
751  // CGV cgxx(cggfs,0.0);
752  // CGV cgdd(cggfs,1.0);
753  // Dune::InverseOperatorResult statstat;
754  // testsolver.apply(native(cgxx),native(cgdd),statstat);
755  // /********/
756 
757  // set up solver
758  int verb=verbose;
759  if (gfs.gridView().comm().rank()>0) verb=0;
760  Solver<V> solver(pop,psp,hybridprec,reduction,maxiter,verb);
761 
762  // solve
763  Dune::InverseOperatorResult stat;
764  watch.reset();
765  solver.apply(z,r,stat);
766  double amg_solve_time = watch.elapsed();
767  if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time << " s" << std::endl;
768  res.converged = stat.converged;
769  res.iterations = stat.iterations;
770  res.elapsed = amg_solve_time+amg_setup_time+triple_product_time;
771  res.reduction = stat.reduction;
772  res.conv_rate = stat.conv_rate;
773  }
774 
775 };
776 }
777 }
778 #endif // DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH
GV::IndexSet IndexSet
Definition: ovlp_amg_dg_backend.hh:135
Dune::PDELab::Backend::Vector< DGGFS, typename DGPrec::domain_type::field_type > V
Definition: ovlp_amg_dg_backend.hh:321
static const int dim
Definition: adaptivity.hh:83
GFS::Traits::GridView GV
Definition: ovlp_amg_dg_backend.hh:134
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: ovlp_amg_dg_backend.hh:561
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
Dune::PDELab::Backend::Vector< CGGFS, field_type > CGV
Definition: ovlp_amg_dg_backend.hh:473
Dune::PDELab::GridOperator< CGGFS, GFS, CGTODGLOP, MBE, field_type, field_type, field_type, CC, CC > PGO
Definition: ovlp_amg_dg_backend.hh:480
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: ovlp_amg_dg_backend.hh:213
Dune::PDELab::Backend::Vector< CGGFS, typename CGPrec::range_type::field_type > CGW
Definition: ovlp_amg_dg_backend.hh:326
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:176
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: ovlp_amg_dg_backend.hh:77
size_t size(EntityType &e) const
Definition: ovlp_amg_dg_backend.hh:70
ISTLBackend_OVLP_AMG_4_DG(DGGO &dggo_, const DGCC &dgcc_, CGGFS &cggfs_, const CGCC &cgcc_, const ParameterTree &params)
Definition: ovlp_amg_dg_backend.hh:618
Dune::PDELab::EmptyTransformation CC
Definition: ovlp_amg_dg_backend.hh:478
Definition: ovlp_amg_dg_backend.hh:304
GV::Grid Grid
Definition: ovlp_amg_dg_backend.hh:42
virtual void apply(V &x, const W &b)
Apply the precondioner.
Definition: ovlp_amg_dg_backend.hh:369
DGGO::Traits::TrialGridFunctionSpace GFS
Definition: ovlp_amg_dg_backend.hh:463
std::map< IndexType, IdType > LocalToGlobalMap
Definition: ovlp_amg_dg_backend.hh:148
OvlpDGAMGPrec(const DGGFS &dggfs_, DGMatrix &dgmatrix_, DGPrec &dgprec_, const DGCC &dgcc_, const CGGFS &cggfs_, CGPrec &cgprec_, const CGCC &cgcc_, P &p_, const DGHelper &dghelper_, const Comm &comm_, int n1_, int n2_)
Constructor.
Definition: ovlp_amg_dg_backend.hh:341
std::map< IdType, IndexType > GlobalToLocalMap
Definition: ovlp_amg_dg_backend.hh:51
Dune::Amg::Parameters Parameters
Definition: ovlp_amg_dg_backend.hh:494
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Backend::Native< CGV > CGVector
Definition: ovlp_amg_dg_backend.hh:474
Backend::Native< PMatrix > P
Definition: ovlp_amg_dg_backend.hh:482
GFS::Traits::GridView GV
Definition: ovlp_amg_dg_backend.hh:39
Backend::Native< W > Y
Definition: ovlp_amg_dg_backend.hh:324
void apply(M &A, V &z, V &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlp_amg_dg_backend.hh:666
GV::Grid Grid
Definition: ovlp_amg_dg_backend.hh:137
Definition: solver.hh:42
Definition: parallelhelper.hh:53
MatrixExchangeDataHandle(const GFS &gfs_, M &m_, const LocalToGlobalMap &l2g_, const GlobalToLocalMap &g2l_, M &m2_)
constructor
Definition: ovlp_amg_dg_backend.hh:244
DGGO::Traits::Jacobian M
Definition: ovlp_amg_dg_backend.hh:466
const std::string s
Definition: function.hh:1101
Dune::BlockPreconditioner< CGVector, CGVector, Comm, Smoother > ParSmoother
Definition: ovlp_amg_dg_backend.hh:492
GV::IndexSet IndexSet
Definition: ovlp_amg_dg_backend.hh:40
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: ovlp_amg_dg_backend.hh:188
M::block_type B
Definition: ovlp_amg_dg_backend.hh:142
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Grid::Traits::GlobalIdSet GlobalIdSet
Definition: ovlp_amg_dg_backend.hh:138
Definition: ovlpistlsolverbackend.hh:389
Definition: genericdatahandle.hh:622
std::map< IndexType, IdType > LocalToGlobalMap
Definition: ovlp_amg_dg_backend.hh:50
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
virtual void post(V &x)
Clean up.
Definition: ovlp_amg_dg_backend.hh:431
std::pair< IdType, B > DataType
export type of data for message buffer
Definition: ovlp_amg_dg_backend.hh:145
Definition: ovlp_amg_dg_backend.hh:127
void jacobian(const Domain &x, Jacobian &a) const
Assembler jacobian.
Definition: gridoperator.hh:205
Dune::PDELab::ISTL::BCRSMatrixBackend MBE
Definition: ovlp_amg_dg_backend.hh:477
Dune::Amg::AMG< ParCGOperator, CGVector, ParSmoother, Comm > AMG
Definition: ovlp_amg_dg_backend.hh:493
Dune::PDELab::Backend::Vector< DGGFS, typename DGPrec::range_type::field_type > W
Definition: ovlp_amg_dg_backend.hh:322
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: ovlp_amg_dg_backend.hh:54
const Entity & e
Definition: localfunctionspace.hh:111
Dune::PDELab::Backend::Matrix< MB, Domain, Range, JF > Jacobian
The type of the jacobian.
Definition: gridoperator.hh:75
Grid::Traits::GlobalIdSet GlobalIdSet
Definition: ovlp_amg_dg_backend.hh:43
sparsity pattern generator
Definition: pattern.hh:13
Dune::SeqSSOR< CGMatrix, CGVector, CGVector, 1 > Smoother
Definition: ovlp_amg_dg_backend.hh:491
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: ovlp_amg_dg_backend.hh:555
Dune::OverlappingSchwarzOperator< CGMatrix, CGVector, CGVector, Comm > ParCGOperator
Definition: ovlp_amg_dg_backend.hh:490
void restore_overlap_entries(const GFS &gfs, Dune::BCRSMatrix< Dune::FieldMatrix< T, n, m >, A > &matrix, Dune::BCRSMatrix< Dune::FieldMatrix< T, n, m >, A > &matrix2)
Definition: ovlp_amg_dg_backend.hh:267
bool fixedsize(int dim, int codim) const
returns true if size per entity of given dim and codim is a constant
Definition: ovlp_amg_dg_backend.hh:60
Definition: ovlp_amg_dg_backend.hh:457
std::map< IdType, IndexType > GlobalToLocalMap
Definition: ovlp_amg_dg_backend.hh:149
int DataType
export type of data for message buffer
Definition: ovlp_amg_dg_backend.hh:47
Definition: ovlp_amg_dg_backend.hh:34
ISTLBackend_OVLP_AMG_4_DG(DGGO &dggo_, const DGCC &dgcc_, CGGFS &cggfs_, const CGCC &cgcc_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: ovlp_amg_dg_backend.hh:574
Backend::Native< V > Vector
Definition: ovlp_amg_dg_backend.hh:469
DGGO::Traits::Domain V
Definition: ovlp_amg_dg_backend.hh:467
Backend::Native< M > Matrix
Definition: ovlp_amg_dg_backend.hh:468
Dune::MatMultMatResult< PTADG, P >::type CGMatrix
Definition: ovlp_amg_dg_backend.hh:486
Default flags for all local operators.
Definition: flags.hh:18
Definition: ovlpistlsolverbackend.hh:44
GlobalIdSet::IdType IdType
Definition: ovlp_amg_dg_backend.hh:44
bool fixedsize(int dim, int codim) const
returns true if size per entity of given dim and codim is a constant
Definition: ovlp_amg_dg_backend.hh:158
TransferLOP CGTODGLOP
Definition: ovlp_amg_dg_backend.hh:479
STL namespace.
IndexSet::IndexType IndexType
Definition: ovlp_amg_dg_backend.hh:41
Vector::field_type field_type
Definition: ovlp_amg_dg_backend.hh:470
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: ovlp_amg_dg_backend.hh:94
Dune::PDELab::ISTL::CommSelector< s, Dune::MPIHelper::isFake >::type Comm
Definition: ovlp_amg_dg_backend.hh:489
LocalGlobalMapDataHandle(const GFS &gfs_, LocalToGlobalMap &l2g_, GlobalToLocalMap &g2l_)
constructor
Definition: ovlp_amg_dg_backend.hh:107
Dune::TransposedMatMultMatResult< P, Matrix >::type PTADG
Definition: ovlp_amg_dg_backend.hh:485
Definition: ovlpistlsolverbackend.hh:438
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: ovlp_amg_dg_backend.hh:567
size_t size(EntityType &e) const
Definition: ovlp_amg_dg_backend.hh:168
Backend::Native< V > X
Definition: ovlp_amg_dg_backend.hh:323
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:417
const P & p
Definition: constraints.hh:147
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:798
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:190
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > &>::type native(T &t)
Definition: backend/interface.hh:192
GlobalIdSet::IdType IdType
Definition: ovlp_amg_dg_backend.hh:139
Definition: constraintstransformation.hh:111
PMatrix & prolongation_matrix()
Definition: ovlp_amg_dg_backend.hh:534
virtual void pre(V &x, W &b)
Prepare the preconditioner.
Definition: ovlp_amg_dg_backend.hh:355
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: ovlp_amg_dg_backend.hh:152
Dune::PDELab::Backend::Vector< CGGFS, typename CGPrec::domain_type::field_type > CGV
Definition: ovlp_amg_dg_backend.hh:325
void setParameters(const Parameters &amg_parameters_)
set AMG parameters
Definition: ovlp_amg_dg_backend.hh:543
PGO::Jacobian PMatrix
Definition: ovlp_amg_dg_backend.hh:481
IndexSet::IndexType IndexType
Definition: ovlp_amg_dg_backend.hh:136