dune-pdelab  2.4-dev
ovlp_amg_dg_backend.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_OVLP_AMG_DG_BACKEND_HH
2 #define DUNE_PDELAB_OVLP_AMG_DG_BACKEND_HH
3 
4 #include <dune/istl/matrixmatrix.hh>
5 
6 #include <dune/grid/common/datahandleif.hh>
7 
16 
17 namespace Dune {
18  namespace PDELab {
19  //***********************************************************
20  // A data handle / function to communicate matrix entries
21  // in the overlap. It turned out that it is actually not
22  // required to do that, but anyway it is there now.
23  //***********************************************************
24 
27  template<class GFS>
29  : public Dune::CommDataHandleIF<LocalGlobalMapDataHandle<GFS>,int>
30  {
31  public:
32  // some types from the grid view
33  typedef typename GFS::Traits::GridView GV;
34  typedef typename GV::IndexSet IndexSet;
35  typedef typename IndexSet::IndexType IndexType;
36  typedef typename GV::Grid Grid;
37  typedef typename Grid::Traits::GlobalIdSet GlobalIdSet;
38  typedef typename GlobalIdSet::IdType IdType;
39 
41  typedef int DataType;
42 
43  // map types
44  typedef std::map<IndexType,IdType> LocalToGlobalMap;
45  typedef std::map<IdType,IndexType> GlobalToLocalMap;
46 
48  bool contains (int dim, int codim) const
49  {
50  return (codim==0);
51  }
52 
54  bool fixedsize (int dim, int codim) const
55  {
56  return true;
57  }
58 
63  template<class EntityType>
64  size_t size (EntityType& e) const
65  {
66  return 1;
67  }
68 
70  template<class MessageBuffer, class EntityType>
71  void gather (MessageBuffer& buff, const EntityType& e) const
72  {
73  // fill map
74  IndexType myindex = indexSet.index(e);
75  IdType myid = globalIdSet.id(e);
76  l2g[myindex] = myid;
77  g2l[myid] = myindex;
78  //std::cout << gv.comm().rank() << ": gather local=" << myindex << " global=" << myid << std::endl;
79 
80  buff.write(0); // this is a dummy, we are not interested in the data
81  }
82 
87  template<class MessageBuffer, class EntityType>
88  void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
89  {
90  DataType x;
91  buff.read(x); // this is a dummy, we are not interested in the data
92 
93  IndexType myindex = indexSet.index(e);
94  IdType myid = globalIdSet.id(e);
95  l2g[myindex] = myid;
96  g2l[myid] = myindex;
97  //std::cout << gv.comm().rank() << ": scatter local=" << myindex << " global=" << myid << std::endl;
98  }
99 
101  LocalGlobalMapDataHandle (const GFS& gfs_, LocalToGlobalMap& l2g_, GlobalToLocalMap& g2l_)
102  : gfs(gfs_), gv(gfs.gridView()), indexSet(gv.indexSet()),
103  grid(gv.grid()), globalIdSet(grid.globalIdSet()),
104  l2g(l2g_), g2l(g2l_)
105  {
106  }
107 
108  private:
109  const GFS& gfs;
110  GV gv;
111  const IndexSet& indexSet;
112  const Grid& grid;
113  const GlobalIdSet& globalIdSet;
114  LocalToGlobalMap& l2g;
115  GlobalToLocalMap& g2l;
116  };
117 
118 
119  // A DataHandle class to exchange rows of a sparse matrix
120  template<class GFS, class M> // mapper type and vector type
122  : public Dune::CommDataHandleIF<MatrixExchangeDataHandle<GFS,M>,
123  std::pair<typename GFS::Traits::GridView::Grid::Traits::GlobalIdSet::IdType,
124  typename M::block_type> >
125  {
126  public:
127  // some types from the grid view
128  typedef typename GFS::Traits::GridView GV;
129  typedef typename GV::IndexSet IndexSet;
130  typedef typename IndexSet::IndexType IndexType;
131  typedef typename GV::Grid Grid;
132  typedef typename Grid::Traits::GlobalIdSet GlobalIdSet;
133  typedef typename GlobalIdSet::IdType IdType;
134 
135  // some types from the matrix
136  typedef typename M::block_type B;
137 
139  typedef std::pair<IdType,B> DataType;
140 
141  // map types
142  typedef std::map<IndexType,IdType> LocalToGlobalMap;
143  typedef std::map<IdType,IndexType> GlobalToLocalMap;
144 
146  bool contains (int dim, int codim) const
147  {
148  return (codim==0);
149  }
150 
152  bool fixedsize (int dim, int codim) const
153  {
154  return false;
155  }
156 
161  template<class EntityType>
162  size_t size (EntityType& e) const
163  {
164  IndexType myindex = indexSet.index(e);
165  typename M::row_type myrow = m[myindex];
166  typename M::row_type::iterator endit=myrow.end();
167  size_t count=0;
168  for (typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it)
169  {
170  typename LocalToGlobalMap::const_iterator find=l2g.find(it.index());
171  if (find!=l2g.end())
172  {
173  count++;
174  }
175  }
176  //std::cout << gv.comm().rank() << ": row=" << myindex << " size=" << count << std::endl;
177  return count;
178  }
179 
181  template<class MessageBuffer, class EntityType>
182  void gather (MessageBuffer& buff, const EntityType& e) const
183  {
184  IndexType myindex = indexSet.index(e);
185  //std::cout << gv.comm().rank() << ": begin gather local=" << myindex << std::endl;
186  typename M::row_type myrow = m[myindex];
187  typename M::row_type::iterator endit=myrow.end();
188  size_t count = 0;
189  for (typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it)
190  {
191  typename LocalToGlobalMap::const_iterator find=l2g.find(it.index());
192  if (find!=l2g.end())
193  {
194  buff.write(std::make_pair(find->second,*it));
195  //std::cout << gv.comm().rank() << ": ==> local=" << find->first << " global=" << find->second << std::endl;
196  count++;
197  }
198  }
199  //std::cout << gv.comm().rank() << ": end gather row=" << myindex << " size=" << count << std::endl;
200  }
201 
206  template<class MessageBuffer, class EntityType>
207  void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
208  {
209  IndexType myindex = indexSet.index(e);
210  std::cout << gv.comm().rank() << ": begin scatter local=" << myindex << " size=" << n << std::endl;
211  DataType x;
212  size_t count = 0;
213  for (size_t i=0; i<n; ++i)
214  {
215  buff.read(x);
216  std::cout << gv.comm().rank() << ": --> received global=" << x.first << std::endl;
217  typename GlobalToLocalMap::const_iterator find=g2l.find(x.first);
218  if (find!=g2l.end())
219  {
220  IndexType nbindex=find->second;
221  if (m.exists(myindex,nbindex))
222  {
223  m[myindex][nbindex] = x.second;
224  B block(x.second);
225  block -= m2[myindex][nbindex];
226  std::cout << gv.comm().rank() << ": compare i=" << myindex << " j=" << nbindex
227  << " norm=" << block.infinity_norm() << std::endl;
228 
229  count++;
230  //std::cout << gv.comm().rank() << ": --> local=" << find->first << " global=" << find->second << std::endl;
231  }
232  }
233  }
234  //std::cout << gv.comm().rank() << ": end scatter row=" << myindex << " size=" << count << std::endl;
235  }
236 
238  MatrixExchangeDataHandle (const GFS& gfs_, M& m_, const LocalToGlobalMap& l2g_, const GlobalToLocalMap& g2l_,M& m2_)
239  : gfs(gfs_), m(m_), gv(gfs.gridView()), indexSet(gv.indexSet()),
240  grid(gv.grid()), globalIdSet(grid.globalIdSet()),
241  l2g(l2g_), g2l(g2l_), m2(m2_)
242  {
243  }
244 
245  private:
246  const GFS& gfs;
247  M& m;
248  GV gv;
249  const IndexSet& indexSet;
250  const Grid& grid;
251  const GlobalIdSet& globalIdSet;
252  const LocalToGlobalMap& l2g;
253  const GlobalToLocalMap& g2l;
254  M& m2;
255  };
256 
257 
260  template<class GFS, class T, class A, int n, int m>
261  void restore_overlap_entries (const GFS& gfs, Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix,
262  Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix2)
263  {
264  typedef Dune::FieldMatrix<T,n,m> B;
265  typedef Dune::BCRSMatrix<B,A> M; // m is of type M
266  typedef typename LocalGlobalMapDataHandle<GFS>::LocalToGlobalMap LocalToGlobalMap;
267  typedef typename LocalGlobalMapDataHandle<GFS>::GlobalToLocalMap GlobalToLocalMap;
268 
269  // build up two maps to associate local indices and global ids
270  LocalToGlobalMap l2g;
271  GlobalToLocalMap g2l;
272  LocalGlobalMapDataHandle<GFS> lgdh(gfs,l2g,g2l);
273  if (gfs.gridView().comm().size()>1)
274  gfs.gridView().communicate(lgdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
275 
276  // exchange matrix entries
277  MatrixExchangeDataHandle<GFS,M> mexdh(gfs,matrix,l2g,g2l,matrix2);
278  if (gfs.gridView().comm().size()>1)
279  gfs.gridView().communicate(mexdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
280  }
281 
282  //***********************************************************
283  // The DG/AMG preconditioner in the overlapping case
284  //***********************************************************
285 
295  template<class DGGFS, class DGMatrix, class DGPrec, class DGCC,
296  class CGGFS, class CGPrec, class CGCC,
297  class P, class DGHelper, class Comm>
299  : public Dune::Preconditioner<Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::domain_type::field_type>,
300  Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::range_type::field_type>>
301  {
302  const DGGFS& dggfs;
303  DGMatrix& dgmatrix;
304  DGPrec& dgprec;
305  const DGCC& dgcc;
306  const CGGFS& cggfs;
307  CGPrec& cgprec;
308  const CGCC& cgcc;
309  P& p;
310  const DGHelper& dghelper;
311  const Comm& comm;
312  int n1,n2;
313 
314  public:
317  typedef typename V::BaseT X;
318  typedef typename W::BaseT Y;
321 
322  // define the category
323  enum {
325  category=Dune::SolverCategory::overlapping
326  };
327 
335  OvlpDGAMGPrec (const DGGFS& dggfs_, DGMatrix& dgmatrix_, DGPrec& dgprec_, const DGCC& dgcc_,
336  const CGGFS& cggfs_, CGPrec& cgprec_, const CGCC& cgcc_, P& p_,
337  const DGHelper& dghelper_, const Comm& comm_, int n1_, int n2_)
338  : dggfs(dggfs_), dgmatrix(dgmatrix_), dgprec(dgprec_), dgcc(dgcc_),
339  cggfs(cggfs_), cgprec(cgprec_), cgcc(cgcc_), p(p_), dghelper(dghelper_),
340  comm(comm_), n1(n1_), n2(n2_)
341  {
342  }
343 
349  virtual void pre (V& x, W& b)
350  {
351  using Backend::native;
352  dgprec.pre(native(x),native(b));
353  CGW cgd(cggfs,0.0);
354  CGV cgv(cggfs,0.0);
355  cgprec.pre(native(cgv),native(cgd));
356  }
357 
363  virtual void apply (V& x, const W& b)
364  {
365  using Backend::native;
366  // need local copies to store defect and solution
367  W d(b);
369  V v(x); // only to get correct size
370 
371  // pre-smoothing on DG matrix
372  for (int i=0; i<n1; i++)
373  {
374  using Backend::native;
375  v = 0.0;
376  dgprec.apply(native(v),native(d));
378  if (dggfs.gridView().comm().size()>1)
379  dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
380  dgmatrix.mmv(native(v),native(d));
382  x += v;
383  }
384 
385  // restrict defect to CG subspace
386  dghelper.maskForeignDOFs(d); // DG defect is additive for overlap 1, but in case we use more
387  CGW cgd(cggfs,0.0);
388  p.mtv(native(d),native(cgd));
390  if (cggfs.gridView().comm().size()>1)
391  cggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication); // now we have consistent defect on coarse grid
393  comm.project(native(cgd));
394  CGV cgv(cggfs,0.0);
395 
396 
397  // call preconditioner
398  cgprec.apply(native(cgv),native(cgd));
399 
400  // prolongate correction
401  p.mv(native(cgv),native(v));
402  dgmatrix.mmv(native(v),native(d));
404  x += v;
405 
406  // post-smoothing on DG matrix
407  for (int i=0; i<n2; i++)
408  {
409  v = 0.0;
410  dgprec.apply(native(v),native(d));
412  if (dggfs.gridView().comm().size()>1)
413  dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
414  dgmatrix.mmv(native(v),native(d));
416  x += v;
417  }
418  }
419 
425  virtual void post (V& x)
426  {
427  dgprec.post(Backend::native(x));
428  CGV cgv(cggfs,0.0);
429  cgprec.post(Backend::native(cgv));
430  }
431  };
432 
433 //***********************************************************
434 // The DG/AMG linear solver backend in the overlapping case
435 //***********************************************************
436 
449 template<class DGGO, class DGCC, class CGGFS, class CGCC, class TransferLOP,
450  template<class,class,class,int> class DGPrec, template<class> class Solver, int s=96>
452  public Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>,
454 {
455 public:
456  // DG grid function space
457  typedef typename DGGO::Traits::TrialGridFunctionSpace GFS;
458 
459  // vectors and matrices on DG level
460  typedef typename DGGO::Traits::Jacobian M; // wrapped istl DG matrix
461  typedef typename DGGO::Traits::Domain V; // wrapped istl DG vector
462  typedef typename M::BaseT Matrix; // istl DG matrix
463  typedef typename V::BaseT Vector; // istl DG vector
464  typedef typename Vector::field_type field_type;
465 
466  // vectors and matrices on CG level
467  using CGV = Dune::PDELab::Backend::Vector<CGGFS,field_type>; // wrapped istl CG vector
468  typedef typename CGV::BaseT CGVector; // istl CG vector
469 
470  // prolongation matrix
473  typedef TransferLOP CGTODGLOP; // local operator
475  typedef typename PGO::Jacobian PMatrix; // wrapped ISTL prolongation matrix
476  typedef typename PMatrix::BaseT P; // ISTL prolongation matrix
477 
478  // CG subspace matrix
479  typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG;
480  typedef typename Dune::MatMultMatResult<PTADG,P>::type CGMatrix; // istl coarse space matrix
481 
482 private:
483 
484  const GFS& gfs;
485  DGGO& dggo;
486  const DGCC& dgcc;
487  CGGFS& cggfs;
488  const CGCC& cgcc;
489  unsigned maxiter;
490  int verbose;
491  bool usesuperlu;
492  std::size_t low_order_space_entries_per_row;
493 
494  CGTODGLOP cgtodglop; // local operator to assemble prolongation matrix
495  PGO pgo; // grid operator to assemble prolongation matrix
496  PMatrix pmatrix; // wrapped prolongation matrix
497 
500  class EmptyLop : public Dune::PDELab::NumericalJacobianApplyVolume<EmptyLop>,
504  {
505  };
506 
507 public:
508 
509  // access to prolongation matrix
511  {
512  return pmatrix;
513  }
514 
517  ISTLBackend_OVLP_AMG_4_DG(DGGO& dggo_, const DGCC& dgcc_, CGGFS& cggfs_, const CGCC& cgcc_,
518  unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
519  : Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>(dggo_.trialGridFunctionSpace())
520  , gfs(dggo_.trialGridFunctionSpace())
521  , dggo(dggo_)
522  , dgcc(dgcc_)
523  , cggfs(cggfs_)
524  , cgcc(cgcc_)
525  , maxiter(maxiter_)
526  , verbose(verbose_)
527  , usesuperlu(usesuperlu_)
528  , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
529  , cgtodglop()
530  , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
531  , pmatrix(pgo)
532  {
533 #if !HAVE_SUPERLU
534  if (usesuperlu == true)
535  {
536  if (gfs.gridView().comm().rank()==0)
537  std::cout << "WARNING: You are using AMG without SuperLU!"
538  << " Please consider installing SuperLU,"
539  << " or set the usesuperlu flag to false"
540  << " to suppress this warning." << std::endl;
541  }
542 #endif
543 
544  // assemble prolongation matrix; this will not change from one apply to the next
545  pmatrix = 0.0;
546  if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
547  CGV cgx(cggfs,0.0); // need vector to call jacobian
548  pgo.jacobian(cgx,pmatrix);
549  }
550 
551 
554  ISTLBackend_OVLP_AMG_4_DG(DGGO& dggo_, const DGCC& dgcc_, CGGFS& cggfs_, const CGCC& cgcc_,
555  const ParameterTree& params)
556  : Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>(dggo_.trialGridFunctionSpace())
557  , gfs(dggo_.trialGridFunctionSpace())
558  , dggo(dggo_)
559  , dgcc(dgcc_)
560  , cggfs(cggfs_)
561  , cgcc(cgcc_)
562  , maxiter(params.get<int>("max_iterations",5000))
563  , verbose(params.get<int>("verbose",1))
564  , usesuperlu(params.get<bool>("use_superlu",true))
565  , low_order_space_entries_per_row(params.get<std::size_t>("low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power))
566  , cgtodglop()
567  , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
568  , pmatrix(pgo)
569  {
570 #if !HAVE_SUPERLU
571  if (usesuperlu == true)
572  {
573  if (gfs.gridView().comm().rank()==0)
574  std::cout << "WARNING: You are using AMG without SuperLU!"
575  << " Please consider installing SuperLU,"
576  << " or set the usesuperlu flag to false"
577  << " to suppress this warning." << std::endl;
578  }
579 #endif
580 
581  // assemble prolongation matrix; this will not change from one apply to the next
582  pmatrix = 0.0;
583  if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
584  CGV cgx(cggfs,0.0); // need vector to call jacobian
585  pgo.jacobian(cgx,pmatrix);
586  }
587 
588 
596  void apply (M& A, V& z, V& r, typename V::ElementType reduction)
597  {
598  using Backend::native;
599  // make operator and scalar product for overlapping solver
601  POP pop(dgcc,A);
603  PSP psp(*this);
604 
605  // compute CG matrix
606  // make grid operator with empty local operator => matrix data type and constraints assembly
607  EmptyLop emptylop;
609  CGGO cggo(cggfs,cgcc,cggfs,cgcc,emptylop,MBE(low_order_space_entries_per_row));
610  typedef typename CGGO::Jacobian CGM;
611 
612  // do triple matrix product ACG = P^T ADG P; this is purely local
613  Dune::Timer watch;
614  watch.reset();
616  CGM acg(attached_container);
617  {
618  PTADG ptadg;
619  Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A)); // 1a
620  //Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A2)); // 1b
621  Dune::matMultMat(native(acg),ptadg,native(pmatrix));
622  }
623  double triple_product_time = watch.elapsed();
624  if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "=== triple matrix product " << triple_product_time << " s" << std::endl;
625  //Dune::printmatrix(std::cout,native(acg),"triple product matrix","row",10,2);
626  CGV cgx(cggfs,0.0); // need vector to call jacobian
627  cggo.jacobian(cgx,acg); // insert trivial rows at processor boundaries
628  //std::cout << "CG constraints: " << cgcc.size() << " out of " << cggfs.globalSize() << std::endl;
629 
630  // NOW we need to insert the processor boundary conditions in DG matrix
632  DGGOEmpty dggoempty(gfs,dgcc,gfs,dgcc,emptylop,MBE(1 << GFS::Traits::GridView::dimension));
633  dggoempty.jacobian(z,A);
634 
635  // and in the residual
637 
638  // now set up parallel AMG solver for the CG subspace
640  Comm oocc(gfs.gridView().comm());
642  CGHELPER cghelper(cggfs,2);
643  cghelper.createIndexSetAndProjectForAMG(acg,oocc);
644  typedef Dune::OverlappingSchwarzOperator<CGMatrix,CGVector,CGVector,Comm> ParCGOperator;
645  ParCGOperator paroop(native(acg),oocc);
646  Dune::OverlappingSchwarzScalarProduct<CGVector,Comm> sp(oocc);
647  typedef Dune::Amg::Parameters Parameters; // AMG parameters (might be nice to change from outside)
648  Parameters params(15,2000);
649  params.setDefaultValuesIsotropic(CGGFS::Traits::GridViewType::Traits::Grid::dimension);
650  params.setDebugLevel(verbose);
651  params.setCoarsenTarget(2000);
652  params.setMaxLevel(20);
653  params.setProlongationDampingFactor(1.6);
654  params.setNoPreSmoothSteps(3);
655  params.setNoPostSmoothSteps(3);
656  params.setGamma(1);
657  params.setAdditive(false);
658  //params.setAccumulate(Dune::Amg::AccumulationMode::noAccu); // atOnceAccu results in deadlock
659  typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother;
660  typedef Dune::BlockPreconditioner<CGVector,CGVector,Comm,Smoother> ParSmoother;
661  typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
662  SmootherArgs smootherArgs;
663  smootherArgs.iterations = 2;
664  smootherArgs.relaxationFactor = 0.92;
665  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
666  Criterion criterion(params);
667  typedef Dune::Amg::AMG<ParCGOperator,CGVector,ParSmoother,Comm> AMG;
668  watch.reset();
669  AMG amg(paroop,criterion,smootherArgs,oocc);
670  double amg_setup_time = watch.elapsed();
671  if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "=== AMG setup " <<amg_setup_time << " s" << std::endl;
672 
673  // set up hybrid DG/CG preconditioner
674  typedef DGPrec<Matrix,Vector,Vector,1> DGPrecType;
675  DGPrecType dgprec(native(A),1,0.92);
676  //DGPrecType dgprec(native(A),0.92);
679  HybridPrec hybridprec(gfs,native(A),dgprec,dgcc,cggfs,amg,cgcc,native(pmatrix),
680  this->parallelHelper(),oocc,3,3);
681 
682  // /********/
683  // /* Test */
684  // /********/
685  // Solver<CGVector> testsolver(paroop,sp,amg,1e-8,100,2);
686  // CGV cgxx(cggfs,0.0);
687  // CGV cgdd(cggfs,1.0);
688  // Dune::InverseOperatorResult statstat;
689  // testsolver.apply(native(cgxx),native(cgdd),statstat);
690  // /********/
691 
692  // set up solver
693  int verb=verbose;
694  if (gfs.gridView().comm().rank()>0) verb=0;
695  Solver<V> solver(pop,psp,hybridprec,reduction,maxiter,verb);
696 
697  // solve
698  Dune::InverseOperatorResult stat;
699  watch.reset();
700  solver.apply(z,r,stat);
701  double amg_solve_time = watch.elapsed();
702  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;
703  res.converged = stat.converged;
704  res.iterations = stat.iterations;
705  res.elapsed = amg_solve_time+amg_setup_time+triple_product_time;
706  res.reduction = stat.reduction;
707  res.conv_rate = stat.conv_rate;
708  }
709 
710 };
711 }
712 }
713 #endif
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: ovlp_amg_dg_backend.hh:146
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:152
GlobalIdSet::IdType IdType
Definition: ovlp_amg_dg_backend.hh:38
void jacobian(const Domain &x, Jacobian &a) const
Assembler jacobian.
Definition: gridoperator.hh:172
Definition: istl/ovlpistlsolverbackend.hh:323
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: ovlp_amg_dg_backend.hh:88
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: ovlp_amg_dg_backend.hh:182
Definition: solver.hh:42
GFS::Traits::GridView GV
Definition: ovlp_amg_dg_backend.hh:128
Dune::PDELab::Backend::Vector< CGGFS, typename CGPrec::range_type::field_type > CGW
Definition: ovlp_amg_dg_backend.hh:320
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:803
IndexSet::IndexType IndexType
Definition: ovlp_amg_dg_backend.hh:130
DGGO::Traits::TrialGridFunctionSpace GFS
Definition: ovlp_amg_dg_backend.hh:457
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: ovlp_amg_dg_backend.hh:48
size_t size(EntityType &e) const
Definition: ovlp_amg_dg_backend.hh:64
virtual void post(V &x)
Clean up.
Definition: ovlp_amg_dg_backend.hh:425
Vector::field_type field_type
Definition: ovlp_amg_dg_backend.hh:464
MatrixExchangeDataHandle(const GFS &gfs_, M &m_, const LocalToGlobalMap &l2g_, const GlobalToLocalMap &g2l_, M &m2_)
constructor
Definition: ovlp_amg_dg_backend.hh:238
GV::IndexSet IndexSet
Definition: ovlp_amg_dg_backend.hh:34
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:335
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:112
const E & e
Definition: interpolate.hh:172
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
RFType reduction
Definition: solver.hh:35
STL namespace.
PGO::Jacobian PMatrix
Definition: ovlp_amg_dg_backend.hh:475
PMatrix::BaseT P
Definition: ovlp_amg_dg_backend.hh:476
Definition: parallelhelper.hh:45
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:261
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: ovlp_amg_dg_backend.hh:71
ISTLBackend_OVLP_AMG_4_DG(DGGO &dggo_, const DGCC &dgcc_, CGGFS &cggfs_, const CGCC &cgcc_, const ParameterTree &params)
Definition: ovlp_amg_dg_backend.hh:554
virtual void pre(V &x, W &b)
Prepare the preconditioner.
Definition: ovlp_amg_dg_backend.hh:349
GV::IndexSet IndexSet
Definition: ovlp_amg_dg_backend.hh:129
std::map< IdType, IndexType > GlobalToLocalMap
Definition: ovlp_amg_dg_backend.hh:45
GlobalIdSet::IdType IdType
Definition: ovlp_amg_dg_backend.hh:133
unsigned int iterations
Definition: solver.hh:33
Backend::attached_container attached_container
Definition: backend/common/tags.hh:37
The category the preconditioner is part of.
Definition: ovlp_amg_dg_backend.hh:325
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
int DataType
export type of data for message buffer
Definition: ovlp_amg_dg_backend.hh:41
Definition: istl/ovlpistlsolverbackend.hh:371
ISTLBackend_OVLP_AMG_4_DG(DGGO &dggo_, const DGCC &dgcc_, CGGFS &cggfs_, const CGCC &cgcc_, unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
Definition: ovlp_amg_dg_backend.hh:517
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:198
TransferLOP CGTODGLOP
Definition: ovlp_amg_dg_backend.hh:473
bool converged
Definition: solver.hh:32
Tag for requesting a vector or matrix container with a pre-attached underlying object.
Definition: backend/common/tags.hh:27
CGV::BaseT CGVector
Definition: ovlp_amg_dg_backend.hh:468
GV::Grid Grid
Definition: ovlp_amg_dg_backend.hh:36
GFS::Traits::GridView GV
Definition: ovlp_amg_dg_backend.hh:33
size_t size(EntityType &e) const
Definition: ovlp_amg_dg_backend.hh:162
Dune::PDELab::Backend::Matrix< MBE, Domain, Range, field_type > Jacobian
The type of the jacobian.
Definition: gridoperator.hh:46
RFType conv_rate
Definition: solver.hh:36
double elapsed
Definition: solver.hh:34
Dune::MatMultMatResult< PTADG, P >::type CGMatrix
Definition: ovlp_amg_dg_backend.hh:480
Definition: genericdatahandle.hh:622
Grid::Traits::GlobalIdSet GlobalIdSet
Definition: ovlp_amg_dg_backend.hh:37
LocalGlobalMapDataHandle(const GFS &gfs_, LocalToGlobalMap &l2g_, GlobalToLocalMap &g2l_)
constructor
Definition: ovlp_amg_dg_backend.hh:101
Grid::Traits::GlobalIdSet GlobalIdSet
Definition: ovlp_amg_dg_backend.hh:132
W::BaseT Y
Definition: ovlp_amg_dg_backend.hh:318
Definition: adaptivity.hh:27
V::BaseT X
Definition: ovlp_amg_dg_backend.hh:317
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:421
Definition: constraintstransformation.hh:111
Dune::PDELab::GridOperator< CGGFS, GFS, CGTODGLOP, MBE, field_type, field_type, field_type, CC, CC > PGO
Definition: ovlp_amg_dg_backend.hh:474
GV::Grid Grid
Definition: ovlp_amg_dg_backend.hh:131
DGGO::Traits::Jacobian M
Definition: ovlp_amg_dg_backend.hh:460
const istl::ParallelHelper< DGGO::Traits::TrialGridFunctionSpace > & parallelHelper() const
Definition: istl/ovlpistlsolverbackend.hh:353
M::block_type B
Definition: ovlp_amg_dg_backend.hh:136
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:54
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:187
static const int dim
Definition: adaptivity.hh:83
Dune::TransposedMatMultMatResult< P, Matrix >::type PTADG
Definition: ovlp_amg_dg_backend.hh:479
Dune::PDELab::Backend::Vector< CGGFS, field_type > CGV
Definition: ovlp_amg_dg_backend.hh:467
Definition: istl/ovlpistlsolverbackend.hh:41
M::BaseT Matrix
Definition: ovlp_amg_dg_backend.hh:462
const std::string s
Definition: function.hh:1102
Definition: ovlp_amg_dg_backend.hh:28
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: ovlp_amg_dg_backend.hh:207
Definition: ovlp_amg_dg_backend.hh:298
V::BaseT Vector
Definition: ovlp_amg_dg_backend.hh:463
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:52
PMatrix & prolongation_matrix()
Definition: ovlp_amg_dg_backend.hh:510
Dune::PDELab::Backend::Vector< CGGFS, typename CGPrec::domain_type::field_type > CGV
Definition: ovlp_amg_dg_backend.hh:319
Dune::PDELab::Backend::Vector< DGGFS, typename DGPrec::domain_type::field_type > V
Definition: ovlp_amg_dg_backend.hh:315
std::pair< IdType, B > DataType
export type of data for message buffer
Definition: ovlp_amg_dg_backend.hh:139
Dune::PDELab::istl::BCRSMatrixBackend MBE
Definition: ovlp_amg_dg_backend.hh:471
Definition: ovlp_amg_dg_backend.hh:121
Dune::PDELab::EmptyTransformation CC
Definition: ovlp_amg_dg_backend.hh:472
std::map< IndexType, IdType > LocalToGlobalMap
Definition: ovlp_amg_dg_backend.hh:142
DGGO::Traits::Domain V
Definition: ovlp_amg_dg_backend.hh:461
virtual void apply(V &x, const W &b)
Apply the precondioner.
Definition: ovlp_amg_dg_backend.hh:363
Dune::PDELab::Backend::Vector< DGGFS, typename DGPrec::range_type::field_type > W
Definition: ovlp_amg_dg_backend.hh:316
std::map< IndexType, IdType > LocalToGlobalMap
Definition: ovlp_amg_dg_backend.hh:44
IndexSet::IndexType IndexType
Definition: ovlp_amg_dg_backend.hh:35
Definition: ovlp_amg_dg_backend.hh:451
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
std::map< IdType, IndexType > GlobalToLocalMap
Definition: ovlp_amg_dg_backend.hh:143
void apply(M &A, V &z, V &r, typename V::ElementType reduction)
solve the given linear system
Definition: ovlp_amg_dg_backend.hh:596